Thermal Tolerance with Color Score and Bleaching Metrics

Setup

Load Packages

##Install Packages if Needed
if (!require("ggplot2")) install.packages("ggplot2")

Loading required package: ggplot2

if (!require("effectsize")) install.packages("effectsize")

Loading required package: effectsize

if (!require("emmeans")) install.packages("emmeans")

Loading required package: emmeans

if (!require("dplyr")) install.packages("dplyr")

Loading required package: dplyr

Attaching package: ‘dplyr’

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

if (!require("tidyr")) install.packages("tidyr")

Loading required package: tidyr

if (!require("Rmisc")) install.packages("Rmisc")

Loading required package: Rmisc
Loading required package: lattice
Loading required package: plyr
----------------------------------------------------------------------------------
You have loaded plyr after dplyr - this is likely to cause problems.
If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
library(plyr); library(dplyr)
----------------------------------------------------------------------------------

Attaching package: ‘plyr’

The following objects are masked from ‘package:dplyr’:

    arrange, count, desc, failwith, id, mutate, rename, summarise,
    summarize

if (!require("ggpubr")) install.packages("ggpubr")

Loading required package: ggpubr

Attaching package: ‘ggpubr’

The following object is masked from ‘package:plyr’:

    mutate

if (!require("cowplot")) install.packages("cowplot")

Loading required package: cowplot

Attaching package: ‘cowplot’

The following object is masked from ‘package:ggpubr’:

    get_legend

if (!require("gridGraphics")) install.packages("gridGraphics")

Loading required package: gridGraphics
Warning: package ‘gridGraphics’ was built under R version 4.3.2Loading required package: grid

if (!require("tidyr")) install.packages("tidyr")

##Load Packages
library(ggplot2) #Required for plotting
library(effectsize) #Required for eta_squared effect sizes
library(emmeans) #Required for pairwise comparisons
library(dplyr) #Required to seperate columns in dataframe
library(tidyr) #Required for data organization
library(Rmisc) #Required for summarySE for summary statistics
library(ggpubr) #Required for adding pairwise p-values to plots with stat_pvalue_manual 
library(cowplot) #Required for arranging ggplots 
library(gridGraphics) #Required for adding labels to arranged plots
library(tidyr) #Required for reshaping datafrom from a wide to long format.

Note: Run “Graphing Parameters” section from 01_ExperimentalSetup.R file

Load and Organize Data

Note: Full Data with Bleaching Metrics and Color Scores created in 04_Models.R file

#Load Data
FullData<-read.csv("Outputs/FullData.csv", header=TRUE)


#Set factor variables 
FullData$TimeP<-factor(FullData$TimeP, levels=c("W2", "M1", "M4"), ordered=TRUE)
FullData$Site<-factor(FullData$Site, levels=c("KL", "SS"), ordered=TRUE)
FullData$Genotype<-factor(FullData$Genotype, levels=c("AC10", "AC12", "AC8"), ordered=TRUE)
FullData$Treatment<-factor(FullData$Treatment, levels=c("Control", "Heat"), ordered=TRUE)
FullData$Treat<-factor(FullData$Treat, levels=c("C", "H"), ordered=TRUE)
FullData$Seas<-factor(FullData$Seas, levels=c("Summer1", "Summer2", "Winter"), ordered=TRUE)

Percent Retention

Subset Thermal Assay Data by Treatment

##Control
FullData_C<-subset(FullData, Treat=="C")

##Heated
FullData_H<-subset(FullData, Treat=="H")

Calculate Percent Retention

Calculating retention as proportion remaining relative to corresponding control levels (0-1) for each site and genotype at each timepoint.

##Calculate averages for Control Treatment for each Site, Genotype, Timepoint, and Season
names(FullData_C)

 [1] "ID"          "RandN"       "TimeP"       "Site"        "Genotype"   
 [6] "Treat"       "Treatment"   "Seas"        "Set"         "Score_Full" 
[11] "Score_TP"    "SA_cm2"      "Chl_ug.cm2"  "Sym10.6_cm2"

FullData_C.a<-aggregate(FullData_C[,c(10:11, 13:14)], list(FullData_C$Site, FullData_C$Genotype, FullData_C$TimeP, FullData_C$Seas), mean, na.action = na.omit)
names(FullData_C.a)[1:4]<-c("Site", "Genotype", "TimeP", "Seas")
names(FullData_C.a)[5:8]<-paste(names(FullData_C)[c(10:11, 13:14)], "C", sep="_")

##Merge Control Averages with Heated Samples
#Merges by Site, Genotype, Timepoint, and Season
TolData<-merge(FullData_H, FullData_C.a, all.x=TRUE )

##Calculate Proportion Retained for each Bleaching Metric relative to Control
TolData$Score_Full.prop<-round((TolData$Score_Full/TolData$Score_Full_C), 4)
TolData$Score_TP.prop<-round((TolData$Score_TP/TolData$Score_TP_C), 4)
TolData$Chl.prop<-round((TolData$Chl_ug.cm2/TolData$Chl_ug.cm2_C), 4)
TolData$Sym.prop<-round((TolData$Sym10.6_cm2/TolData$Sym10.6_cm2_C), 4)

##Set values >1 to 1
TolData$Score_Full.prop[which(TolData$Score_Full.prop>1)]<-1.0000
TolData$Score_TP.prop[which(TolData$Score_TP.prop>1)]<-1.0000
TolData$Chl.prop[which(TolData$Chl.prop>1)]<-1.0000
TolData$Sym.prop[which(TolData$Sym.prop>1)]<-1.0000

##Create Site and Genotype Variable
TolData$Site.Geno<-paste(TolData$Site, TolData$Genotype, sep="_")

Subset by Timepoint

TolData_W2<-subset(TolData, TimeP=="W2")
TolData_M1<-subset(TolData, TimeP=="M1")
TolData_M4<-subset(TolData, TimeP=="M4")

Thermal Tolerance Chlorophyll

Full Set

Run Model

##Check normality
hist(TolData$Chl.prop)

shapiro.test(TolData$Chl.prop)

    Shapiro-Wilk normality test

data:  TolData$Chl.prop
W = 0.97395, p-value = 0.165

#Normal

##Model as a function of Site and Genotype and Timepoint
Tol_Chl_lm<-lm(Chl.prop~Site+Genotype+TimeP+ Site:Genotype + Site:TimeP + Genotype:TimeP, data=TolData)

Check Residuals

##Check Normality of Residuals
#Distribution 
plot(density(resid(Tol_Chl_lm)))

#Q-Q plot
qqnorm(resid(Tol_Chl_lm)); qqline(resid(Tol_Chl_lm))

##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_Chl_lm), resid(Tol_Chl_lm))

Model Results

#Model Results
summary(Tol_Chl_lm)

Call:
lm(formula = Chl.prop ~ Site + Genotype + TimeP + Site:Genotype + 
    Site:TimeP + Genotype:TimeP, data = TolData)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.127060 -0.031006 -0.004219  0.036829  0.097690 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)         0.2100753  0.0065828  31.913  < 2e-16 ***
Site.L             -0.0093229  0.0092942  -1.003 0.320289    
Genotype.L         -0.0980285  0.0111500  -8.792 5.31e-12 ***
Genotype.Q         -0.0467859  0.0116301  -4.023 0.000180 ***
TimeP.L             0.0638187  0.0112994   5.648 6.22e-07 ***
TimeP.Q             0.0464233  0.0115033   4.036 0.000173 ***
Site.L:Genotype.L   0.0003917  0.0157685   0.025 0.980272    
Site.L:Genotype.Q   0.0592529  0.0164183   3.609 0.000673 ***
Site.L:TimeP.L     -0.0528050  0.0159798  -3.304 0.001694 ** 
Site.L:TimeP.Q     -0.0236749  0.0161894  -1.462 0.149435    
Genotype.L:TimeP.L  0.0406686  0.0194339   2.093 0.041093 *  
Genotype.Q:TimeP.L -0.0627603  0.0196595  -3.192 0.002354 ** 
Genotype.L:TimeP.Q  0.0800244  0.0191901   4.170 0.000111 ***
Genotype.Q:TimeP.Q -0.1328570  0.0206170  -6.444 3.28e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.05393 on 54 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.8229,    Adjusted R-squared:  0.7803 
F-statistic:  19.3 on 13 and 54 DF,  p-value: 8.612e-16

anova(Tol_Chl_lm)

Analysis of Variance Table

Response: Chl.prop
               Df   Sum Sq  Mean Sq F value    Pr(>F)    
Site            1 0.001563 0.001563  0.5375  0.466633    
Genotype        2 0.297307 0.148654 51.1117 3.495e-13 ***
TimeP           2 0.132285 0.066142 22.7418 6.845e-08 ***
Site:Genotype   2 0.041114 0.020557  7.0681  0.001877 ** 
Site:TimeP      2 0.038643 0.019322  6.6434  0.002634 ** 
Genotype:TimeP  4 0.218951 0.054738 18.8205 9.691e-10 ***
Residuals      54 0.157054 0.002908                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#Effect Size of Predictors
eta_squared(Tol_Chl_lm, partial=FALSE)

# Effect Size for ANOVA (Type I)

Parameter      |     Eta2 |       95% CI
----------------------------------------
Site           | 1.76e-03 | [0.00, 1.00]
Genotype       |     0.34 | [0.16, 1.00]
TimeP          |     0.15 | [0.02, 1.00]
Site:Genotype  |     0.05 | [0.00, 1.00]
Site:TimeP     |     0.04 | [0.00, 1.00]
Genotype:TimeP |     0.25 | [0.06, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
Tol_Chl_lm.res<-data.frame(anova(Tol_Chl_lm))
Tol_Chl_lm.res$Predictor<-rownames(Tol_Chl_lm.res)
Tol_Chl_lm.res$EtaSq<-c(eta_squared(Tol_Chl_lm, partial=FALSE)$Eta2, NA)
Tol_Chl_lm.res$Response<-rep("Chlorophyll", nrow(Tol_Chl_lm.res))
Tol_Chl_lm.res<-Tol_Chl_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")

Strong influence of Timepoint, including interactions with main variables of interest, Site and Genotype. Will analyze each Timepoint individually.

Chl Summer 1 Timepoint W2

Run Model

##Check normality
hist(TolData_W2$Chl.prop)

shapiro.test(TolData_W2$Chl.prop)

    Shapiro-Wilk normality test

data:  TolData_W2$Chl.prop
W = 0.98407, p-value = 0.9671

#Normal

##Model as a function of Site and Genotype
Tol_Chl_W2_lm<-lm(Chl.prop~Site+Genotype+Site:Genotype, data=TolData_W2)

Check Residuals

##Check Normality of Residuals
#Distribution 
plot(density(resid(Tol_Chl_W2_lm)))

#Q-Q plot
qqnorm(resid(Tol_Chl_W2_lm)); qqline(resid(Tol_Chl_W2_lm))

##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_Chl_W2_lm), resid(Tol_Chl_W2_lm))

Model Results

Overall

#Model Results
summary(Tol_Chl_W2_lm)

Call:
lm(formula = Chl.prop ~ Site + Genotype + Site:Genotype, data = TolData_W2)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.089900 -0.031587 -0.008183  0.038881  0.103550 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)        0.1842264  0.0119450  15.423 5.03e-11 ***
Site.L             0.0188110  0.0168928   1.114 0.281921    
Genotype.L        -0.0947759  0.0204291  -4.639 0.000273 ***
Genotype.Q        -0.0585870  0.0209464  -2.797 0.012921 *  
Site.L:Genotype.L -0.0002833  0.0288911  -0.010 0.992297    
Site.L:Genotype.Q  0.0366569  0.0296227   1.237 0.233769    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.05552 on 16 degrees of freedom
  (1 observation deleted due to missingness)
Multiple R-squared:  0.6758,    Adjusted R-squared:  0.5745 
F-statistic: 6.671 on 5 and 16 DF,  p-value: 0.001551

anova(Tol_Chl_W2_lm)

Analysis of Variance Table

Response: Chl.prop
              Df   Sum Sq  Mean Sq F value   Pr(>F)    
Site           1 0.005580 0.005580  1.8106 0.197195    
Genotype       2 0.092482 0.046241 15.0038 0.000214 ***
Site:Genotype  2 0.004732 0.002366  0.7677 0.480448    
Residuals     16 0.049311 0.003082                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#Effect Size of Predictors
eta_squared(Tol_Chl_W2_lm, partial=FALSE)

# Effect Size for ANOVA (Type I)

Parameter     | Eta2 |       95% CI
-----------------------------------
Site          | 0.04 | [0.00, 1.00]
Genotype      | 0.61 | [0.29, 1.00]
Site:Genotype | 0.03 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
Tol_Chl_W2_lm.res<-data.frame(anova(Tol_Chl_W2_lm))
Tol_Chl_W2_lm.res$Predictor<-rownames(Tol_Chl_W2_lm.res)
Tol_Chl_W2_lm.res$EtaSq<-c(eta_squared(Tol_Chl_W2_lm, partial=FALSE)$Eta2, NA)
Tol_Chl_W2_lm.res$Response<-rep("Chlorophyll", nrow(Tol_Chl_W2_lm.res))
Tol_Chl_W2_lm.res$TimeP<-rep("W2", nrow(Tol_Chl_W2_lm.res))
Tol_Chl_W2_lm.res<-Tol_Chl_W2_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")

Significant effect of Genotype. Still checking Site*Genotype for comparability across Timepoints.

Pairwise

#Pairwise comparisons across:

#Genotypes within Sites
emmeans(Tol_Chl_W2_lm, pairwise~Genotype | Site)

$emmeans
Site = KL:
 Genotype emmean     SE df lower.CL upper.CL
 AC10     0.2033 0.0278 16   0.1445    0.262
 AC12     0.2399 0.0278 16   0.1811    0.299
 AC8      0.0696 0.0278 16   0.0107    0.128

Site = SS:
 Genotype emmean     SE df lower.CL upper.CL
 AC10     0.2514 0.0278 16   0.1925    0.310
 AC12     0.2242 0.0321 16   0.1563    0.292
 AC8      0.1170 0.0321 16   0.0491    0.185

Confidence level used: 0.95 

$contrasts
Site = KL:
 contrast    estimate     SE df t.ratio p.value
 AC10 - AC12  -0.0366 0.0393 16  -0.933  0.6280
 AC10 - AC8    0.1338 0.0393 16   3.407  0.0095
 AC12 - AC8    0.1704 0.0393 16   4.340  0.0014

Site = SS:
 contrast    estimate     SE df t.ratio p.value
 AC10 - AC12   0.0272 0.0424 16   0.640  0.8003
 AC10 - AC8    0.1343 0.0424 16   3.168  0.0156
 AC12 - AC8    0.1072 0.0453 16   2.364  0.0753

P value adjustment: tukey method for comparing a family of 3 estimates 

#Sites within Genotypes
emmeans(Tol_Chl_W2_lm, pairwise~Site | Genotype)

$emmeans
Genotype = AC10:
 Site emmean     SE df lower.CL upper.CL
 KL   0.2033 0.0278 16   0.1445    0.262
 SS   0.2514 0.0278 16   0.1925    0.310

Genotype = AC12:
 Site emmean     SE df lower.CL upper.CL
 KL   0.2399 0.0278 16   0.1811    0.299
 SS   0.2242 0.0321 16   0.1563    0.292

Genotype = AC8:
 Site emmean     SE df lower.CL upper.CL
 KL   0.0696 0.0278 16   0.0107    0.128
 SS   0.1170 0.0321 16   0.0491    0.185

Confidence level used: 0.95 

$contrasts
Genotype = AC10:
 contrast estimate     SE df t.ratio p.value
 KL - SS   -0.0481 0.0393 16  -1.224  0.2387

Genotype = AC12:
 contrast estimate     SE df t.ratio p.value
 KL - SS    0.0157 0.0424 16   0.371  0.7156

Genotype = AC8:
 contrast estimate     SE df t.ratio p.value
 KL - SS   -0.0475 0.0424 16  -1.120  0.2793

##Save p-values

#Genotypes within Sites
Tol_Chl_W2_lm.geno<-data.frame(emmeans(Tol_Chl_W2_lm, pairwise~Genotype | Site)$contrasts)
Tol_Chl_W2_lm.geno<-Tol_Chl_W2_lm.geno %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_Chl_W2_lm.geno$group1<-paste(Tol_Chl_W2_lm.geno$Site, Tol_Chl_W2_lm.geno$group1, sep="_")
Tol_Chl_W2_lm.geno$group2<-paste(Tol_Chl_W2_lm.geno$Site, Tol_Chl_W2_lm.geno$group2, sep="_")

#Sites within Genotypes
Tol_Chl_W2_lm.site<-data.frame(emmeans(Tol_Chl_W2_lm, pairwise~Site | Genotype)$contrasts)
Tol_Chl_W2_lm.site<-Tol_Chl_W2_lm.site %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_Chl_W2_lm.site$group1<-paste(Tol_Chl_W2_lm.site$group1, Tol_Chl_W2_lm.site$Genotype, sep="_")
Tol_Chl_W2_lm.site$group2<-paste(Tol_Chl_W2_lm.site$group2, Tol_Chl_W2_lm.site$Genotype, sep="_")

#Full list of p-values
Tol_Chl_W2_lm.p<-rbind(Tol_Chl_W2_lm.geno[,c(1:2,4:8)], Tol_Chl_W2_lm.site[,c(1:2,4:8)])
Tol_Chl_W2_lm.p<-Tol_Chl_W2_lm.p %>% dplyr::rename( p = p.value)

#Add Significance Levels
Tol_Chl_W2_lm.p$Sig<-ifelse(Tol_Chl_W2_lm.p$p<0.001, "***", ifelse(Tol_Chl_W2_lm.p$p<0.01, "**", ifelse(Tol_Chl_W2_lm.p$p<0.05, "*", NA)))

#Specify Response and Timepoint
Tol_Chl_W2_lm.p$Response<-rep("Chlorophyll", nrow(Tol_Chl_W2_lm.p))
Tol_Chl_W2_lm.p$TimeP<-rep("W2", nrow(Tol_Chl_W2_lm.p))

Plot Retention by Site and Genotype

##Summary statistics by Site and Genotype
Tol_Chl_W2_SG<-summarySE(TolData_W2, measurevar="Chl.prop", groupvars=c("Site.Geno", "Site", "Genotype"), na.rm=TRUE)

##Plot Average Retention across Treatments
Tol_Chl_W2_SG.plot<-ggplot(Tol_Chl_W2_SG, aes(x=Site.Geno, y=Chl.prop, colour=Genotype)) + 
  geom_errorbar(aes(ymin=Chl.prop-se, ymax=Chl.prop+se), width=cap.sz, linewidth=bar.sz)+
  geom_point(size=point.sz)+ 
   ggtitle("Chlorophyll Retention")+
  theme_classic()+
  theme( axis.title.x = element_text(size = axis.title.sz), 
         axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"),
        axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz),
        legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"),
        legend.position="bottom", 
        legend.direction="horizontal",
        plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
  labs(x="", y="Proportion Retained")+
  ylim(0, 1)+ 
  scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
  scale_color_manual(values = Geno.colors.o)+ 
  stat_pvalue_manual(data=Tol_Chl_W2_lm.p,  y.position=0.4, step.increase=0.25, label="Sig", hide.ns=TRUE); Tol_Chl_W2_SG.plot

Chl Summer 2 Timepoint M1

Run Model

##Check normality
hist(TolData_M1$Chl.prop)

shapiro.test(TolData_M1$Chl.prop)

    Shapiro-Wilk normality test

data:  TolData_M1$Chl.prop
W = 0.9306, p-value = 0.1263

#Normal

##Model as a function of Site and Genotype
Tol_Chl_M1_lm<-lm(Chl.prop~Site+Genotype+Site:Genotype, data=TolData_M1)

Check Residuals

##Check Normality of Residuals
#Distribution 
plot(density(resid(Tol_Chl_M1_lm)))

#Q-Q plot
qqnorm(resid(Tol_Chl_M1_lm)); qqline(resid(Tol_Chl_M1_lm))

##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_Chl_M1_lm), resid(Tol_Chl_M1_lm))

Model Results

Overall

#Model Results
summary(Tol_Chl_M1_lm)

Call:
lm(formula = Chl.prop ~ Site + Genotype + Site:Genotype, data = TolData_M1)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.107050 -0.034575  0.003479  0.025531  0.088050 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)        0.17217    0.01178  14.619 1.12e-10 ***
Site.L             0.01036    0.01666   0.622   0.5428    
Genotype.L        -0.16337    0.01935  -8.442 2.74e-07 ***
Genotype.Q         0.06169    0.02139   2.883   0.0108 *  
Site.L:Genotype.L  0.03819    0.02737   1.395   0.1820    
Site.L:Genotype.Q  0.05454    0.03026   1.803   0.0903 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.05474 on 16 degrees of freedom
Multiple R-squared:  0.8424,    Adjusted R-squared:  0.7932 
F-statistic: 17.11 on 5 and 16 DF,  p-value: 6.411e-06

anova(Tol_Chl_M1_lm)

Analysis of Variance Table

Response: Chl.prop
              Df   Sum Sq  Mean Sq F value    Pr(>F)    
Site           1 0.002283 0.002283  0.7619    0.3956    
Genotype       2 0.238424 0.119212 39.7888 6.168e-07 ***
Site:Genotype  2 0.015569 0.007785  2.5982    0.1054    
Residuals     16 0.047938 0.002996                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#Effect Size of Predictors
eta_squared(Tol_Chl_M1_lm, partial=FALSE)

# Effect Size for ANOVA (Type I)

Parameter     |     Eta2 |       95% CI
---------------------------------------
Site          | 7.50e-03 | [0.00, 1.00]
Genotype      |     0.78 | [0.58, 1.00]
Site:Genotype |     0.05 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
Tol_Chl_M1_lm.res<-data.frame(anova(Tol_Chl_M1_lm))
Tol_Chl_M1_lm.res$Predictor<-rownames(Tol_Chl_M1_lm.res)
Tol_Chl_M1_lm.res$EtaSq<-c(eta_squared(Tol_Chl_M1_lm, partial=FALSE)$Eta2, NA)
Tol_Chl_M1_lm.res$Response<-rep("Chlorophyll", nrow(Tol_Chl_M1_lm.res))
Tol_Chl_M1_lm.res$TimeP<-rep("M1", nrow(Tol_Chl_M1_lm.res))
Tol_Chl_M1_lm.res<-Tol_Chl_M1_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")

Significant effect of Genotype. Still checking Site*Genotype for comparability across Timepoints.

Pairwise

#Pairwise comparisons across:

#Genotypes within Sites
emmeans(Tol_Chl_M1_lm, pairwise~Genotype | Site)

$emmeans
Site = KL:
 Genotype emmean     SE df lower.CL upper.CL
 AC10     0.3089 0.0274 16   0.2509   0.3669
 AC12     0.1460 0.0316 16   0.0790   0.2130
 AC8      0.0397 0.0274 16  -0.0183   0.0977

Site = SS:
 Genotype emmean     SE df lower.CL upper.CL
 AC10     0.3169 0.0274 16   0.2588   0.3749
 AC12     0.0976 0.0316 16   0.0306   0.1646
 AC8      0.1240 0.0274 16   0.0660   0.1820

Confidence level used: 0.95 

$contrasts
Site = KL:
 contrast    estimate     SE df t.ratio p.value
 AC10 - AC12   0.1629 0.0418 16   3.897  0.0035
 AC10 - AC8    0.2692 0.0387 16   6.956  <.0001
 AC12 - AC8    0.1063 0.0418 16   2.543  0.0538

Site = SS:
 contrast    estimate     SE df t.ratio p.value
 AC10 - AC12   0.2192 0.0418 16   5.244  0.0002
 AC10 - AC8    0.1928 0.0387 16   4.983  0.0004
 AC12 - AC8   -0.0264 0.0418 16  -0.631  0.8056

P value adjustment: tukey method for comparing a family of 3 estimates 

#Sites within Genotypes
emmeans(Tol_Chl_M1_lm, pairwise~Site | Genotype)

$emmeans
Genotype = AC10:
 Site emmean     SE df lower.CL upper.CL
 KL   0.3089 0.0274 16   0.2509   0.3669
 SS   0.3169 0.0274 16   0.2588   0.3749

Genotype = AC12:
 Site emmean     SE df lower.CL upper.CL
 KL   0.1460 0.0316 16   0.0790   0.2130
 SS   0.0976 0.0316 16   0.0306   0.1646

Genotype = AC8:
 Site emmean     SE df lower.CL upper.CL
 KL   0.0397 0.0274 16  -0.0183   0.0977
 SS   0.1240 0.0274 16   0.0660   0.1820

Confidence level used: 0.95 

$contrasts
Genotype = AC10:
 contrast estimate     SE df t.ratio p.value
 KL - SS  -0.00795 0.0387 16  -0.205  0.8398

Genotype = AC12:
 contrast estimate     SE df t.ratio p.value
 KL - SS   0.04833 0.0447 16   1.081  0.2955

Genotype = AC8:
 contrast estimate     SE df t.ratio p.value
 KL - SS  -0.08432 0.0387 16  -2.179  0.0447

##Save p-values

#Genotypes within Sites
Tol_Chl_M1_lm.geno<-data.frame(emmeans(Tol_Chl_M1_lm, pairwise~Genotype | Site)$contrasts)
Tol_Chl_M1_lm.geno<-Tol_Chl_M1_lm.geno %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_Chl_M1_lm.geno$group1<-paste(Tol_Chl_M1_lm.geno$Site, Tol_Chl_M1_lm.geno$group1, sep="_")
Tol_Chl_M1_lm.geno$group2<-paste(Tol_Chl_M1_lm.geno$Site, Tol_Chl_M1_lm.geno$group2, sep="_")

#Sites within Genotypes
Tol_Chl_M1_lm.site<-data.frame(emmeans(Tol_Chl_M1_lm, pairwise~Site | Genotype)$contrasts)
Tol_Chl_M1_lm.site<-Tol_Chl_M1_lm.site %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_Chl_M1_lm.site$group1<-paste(Tol_Chl_M1_lm.site$group1, Tol_Chl_M1_lm.site$Genotype, sep="_")
Tol_Chl_M1_lm.site$group2<-paste(Tol_Chl_M1_lm.site$group2, Tol_Chl_M1_lm.site$Genotype, sep="_")

#Full list of p-values
Tol_Chl_M1_lm.p<-rbind(Tol_Chl_M1_lm.geno[,c(1:2,4:8)], Tol_Chl_M1_lm.site[,c(1:2,4:8)])
Tol_Chl_M1_lm.p<-Tol_Chl_M1_lm.p %>% dplyr::rename( p = p.value)

#Add Significance Levels
Tol_Chl_M1_lm.p$Sig<-ifelse(Tol_Chl_M1_lm.p$p<0.001, "***", ifelse(Tol_Chl_M1_lm.p$p<0.01, "**", ifelse(Tol_Chl_M1_lm.p$p<0.05, "*", NA)))

#Specify Response and Timepoint
Tol_Chl_M1_lm.p$Response<-rep("Chlorophyll", nrow(Tol_Chl_M1_lm.p))
Tol_Chl_M1_lm.p$TimeP<-rep("M1", nrow(Tol_Chl_M1_lm.p))

Plot Retention by Site and Genotype

##Summary statistics by Site and Genotype
Tol_Chl_M1_SG<-summarySE(TolData_M1, measurevar="Chl.prop", groupvars=c("Site.Geno", "Site", "Genotype"), na.rm=TRUE)

##Plot Average Retention across Treatments
Tol_Chl_M1_SG.plot<-ggplot(Tol_Chl_M1_SG, aes(x=Site.Geno, y=Chl.prop, colour=Genotype)) + 
  geom_errorbar(aes(ymin=Chl.prop-se, ymax=Chl.prop+se), width=cap.sz, linewidth=bar.sz)+
  geom_point(size=point.sz)+ 
   ggtitle("Chlorophyll Retention")+
  theme_classic()+
  theme( axis.title.x = element_text(size = axis.title.sz), 
         axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"),
        axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz),
        legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"),
        legend.position="bottom", 
        legend.direction="horizontal",
        plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
  labs(x="", y="Proportion Retained")+
  ylim(0, 1)+ 
  scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
  scale_color_manual(values = Geno.colors.o)+ 
  stat_pvalue_manual(data=Tol_Chl_M1_lm.p,  y.position=0.4, step.increase=0.25, label="Sig", hide.ns=TRUE); Tol_Chl_M1_SG.plot

Chl Winter Timepoint M4

Run Model

##Check normality
hist(TolData_M4$Chl.prop)

shapiro.test(TolData_M4$Chl.prop)

    Shapiro-Wilk normality test

data:  TolData_M4$Chl.prop
W = 0.91299, p-value = 0.04096

#Not Normal

##Try log+1 transformation
hist(log(TolData_M4$Chl.prop+1))

shapiro.test(log(TolData_M4$Chl.prop+1))

    Shapiro-Wilk normality test

data:  log(TolData_M4$Chl.prop + 1)
W = 0.93102, p-value = 0.1028

#Normal

##Model as a function of Site and Genotype
##Model with log transformation
Tol_Chl_M4_lm<-lm(log(Chl.prop+1)~Site+Genotype+Site:Genotype, data=TolData_M4)

Check Residuals

##Check Normality of Residuals
#Distribution 
plot(density(resid(Tol_Chl_M4_lm)))

#Q-Q plot
qqnorm(resid(Tol_Chl_M4_lm)); qqline(resid(Tol_Chl_M4_lm))

##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_Chl_M4_lm), resid(Tol_Chl_M4_lm))

Model Results

Overall

#Model Results
summary(Tol_Chl_M4_lm)

Call:
lm(formula = log(Chl.prop + 1) ~ Site + Genotype + Site:Genotype, 
    data = TolData_M4)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.062753 -0.025900 -0.004896  0.017009  0.065759 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)        0.23863    0.00807  29.568  < 2e-16 ***
Site.L            -0.04272    0.01141  -3.743  0.00149 ** 
Genotype.L        -0.03027    0.01398  -2.166  0.04400 *  
Genotype.Q        -0.11052    0.01398  -7.907  2.9e-07 ***
Site.L:Genotype.L -0.03189    0.01977  -1.613  0.12415    
Site.L:Genotype.Q  0.05753    0.01977   2.910  0.00934 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.03954 on 18 degrees of freedom
Multiple R-squared:  0.8368,    Adjusted R-squared:  0.7915 
F-statistic: 18.46 on 5 and 18 DF,  p-value: 1.598e-06

anova(Tol_Chl_M4_lm)

Analysis of Variance Table

Response: log(Chl.prop + 1)
              Df   Sum Sq  Mean Sq F value    Pr(>F)    
Site           1 0.021901 0.021901 14.0102  0.001489 ** 
Genotype       2 0.105056 0.052528 33.6031 8.379e-07 ***
Site:Genotype  2 0.017304 0.008652  5.5349  0.013381 *  
Residuals     18 0.028137 0.001563                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#Effect Size of Predictors
eta_squared(Tol_Chl_M4_lm, partial=FALSE)

# Effect Size for ANOVA (Type I)

Parameter     | Eta2 |       95% CI
-----------------------------------
Site          | 0.13 | [0.00, 1.00]
Genotype      | 0.61 | [0.32, 1.00]
Site:Genotype | 0.10 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
Tol_Chl_M4_lm.res<-data.frame(anova(Tol_Chl_M4_lm))
Tol_Chl_M4_lm.res$Predictor<-rownames(Tol_Chl_M4_lm.res)
Tol_Chl_M4_lm.res$EtaSq<-c(eta_squared(Tol_Chl_M4_lm, partial=FALSE)$Eta2, NA)
Tol_Chl_M4_lm.res$Response<-rep("Chlorophyll", nrow(Tol_Chl_M4_lm.res))
Tol_Chl_M4_lm.res$TimeP<-rep("M4", nrow(Tol_Chl_M4_lm.res))
Tol_Chl_M4_lm.res<-Tol_Chl_M4_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")

Significant effect of Site and Genotype and Site*Genotype.

Pairwise

#Pairwise comparisons across:

#Genotypes within Sites
emmeans(Tol_Chl_M4_lm, pairwise~Genotype | Site)

$emmeans
Site = KL:
 Genotype emmean     SE df lower.CL upper.CL
 AC10      0.213 0.0198 18    0.171    0.254
 AC12      0.392 0.0198 18    0.351    0.434
 AC8       0.202 0.0198 18    0.160    0.243

Site = SS:
 Genotype emmean     SE df lower.CL upper.CL
 AC10      0.217 0.0198 18    0.176    0.259
 AC12      0.265 0.0198 18    0.224    0.307
 AC8       0.143 0.0198 18    0.101    0.184

Results are given on the log(mu + 1) (not the response) scale. 
Confidence level used: 0.95 

$contrasts
Site = KL:
 contrast    estimate    SE df t.ratio p.value
 AC10 - AC12  -0.1797 0.028 18  -6.428  <.0001
 AC10 - AC8    0.0109 0.028 18   0.391  0.9196
 AC12 - AC8    0.1906 0.028 18   6.819  <.0001

Site = SS:
 contrast    estimate    SE df t.ratio p.value
 AC10 - AC12  -0.0482 0.028 18  -1.724  0.2237
 AC10 - AC8    0.0747 0.028 18   2.672  0.0393
 AC12 - AC8    0.1229 0.028 18   4.396  0.0010

Note: contrasts are still on the log(mu + 1) scale 
P value adjustment: tukey method for comparing a family of 3 estimates 

#Sites within Genotypes
emmeans(Tol_Chl_M4_lm, pairwise~Site | Genotype)

$emmeans
Genotype = AC10:
 Site emmean     SE df lower.CL upper.CL
 KL    0.213 0.0198 18    0.171    0.254
 SS    0.217 0.0198 18    0.176    0.259

Genotype = AC12:
 Site emmean     SE df lower.CL upper.CL
 KL    0.392 0.0198 18    0.351    0.434
 SS    0.265 0.0198 18    0.224    0.307

Genotype = AC8:
 Site emmean     SE df lower.CL upper.CL
 KL    0.202 0.0198 18    0.160    0.243
 SS    0.143 0.0198 18    0.101    0.184

Results are given on the log(mu + 1) (not the response) scale. 
Confidence level used: 0.95 

$contrasts
Genotype = AC10:
 contrast estimate    SE df t.ratio p.value
 KL - SS  -0.00468 0.028 18  -0.167  0.8688

Genotype = AC12:
 contrast estimate    SE df t.ratio p.value
 KL - SS   0.12684 0.028 18   4.537  0.0003

Genotype = AC8:
 contrast estimate    SE df t.ratio p.value
 KL - SS   0.05909 0.028 18   2.114  0.0488

Note: contrasts are still on the log(mu + 1) scale 

##Save p-values

#Genotypes within Sites
Tol_Chl_M4_lm.geno<-data.frame(emmeans(Tol_Chl_M4_lm, pairwise~Genotype | Site)$contrasts)
Tol_Chl_M4_lm.geno<-Tol_Chl_M4_lm.geno %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_Chl_M4_lm.geno$group1<-paste(Tol_Chl_M4_lm.geno$Site, Tol_Chl_M4_lm.geno$group1, sep="_")
Tol_Chl_M4_lm.geno$group2<-paste(Tol_Chl_M4_lm.geno$Site, Tol_Chl_M4_lm.geno$group2, sep="_")

#Sites within Genotypes
Tol_Chl_M4_lm.site<-data.frame(emmeans(Tol_Chl_M4_lm, pairwise~Site | Genotype)$contrasts)
Tol_Chl_M4_lm.site<-Tol_Chl_M4_lm.site %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_Chl_M4_lm.site$group1<-paste(Tol_Chl_M4_lm.site$group1, Tol_Chl_M4_lm.site$Genotype, sep="_")
Tol_Chl_M4_lm.site$group2<-paste(Tol_Chl_M4_lm.site$group2, Tol_Chl_M4_lm.site$Genotype, sep="_")

#Full list of p-values
Tol_Chl_M4_lm.p<-rbind(Tol_Chl_M4_lm.geno[,c(1:2,4:8)], Tol_Chl_M4_lm.site[,c(1:2,4:8)])
Tol_Chl_M4_lm.p<-Tol_Chl_M4_lm.p %>% dplyr::rename( p = p.value)

#Add Significance Levels
Tol_Chl_M4_lm.p$Sig<-ifelse(Tol_Chl_M4_lm.p$p<0.001, "***", ifelse(Tol_Chl_M4_lm.p$p<0.01, "**", ifelse(Tol_Chl_M4_lm.p$p<0.05, "*", NA)))

#Specify Response and Timepoint
Tol_Chl_M4_lm.p$Response<-rep("Chlorophyll", nrow(Tol_Chl_M4_lm.p))
Tol_Chl_M4_lm.p$TimeP<-rep("M4", nrow(Tol_Chl_M4_lm.p))

Plot Retention by Site and Genotype

##Summary statistics by Site and Genotype
Tol_Chl_M4_SG<-summarySE(TolData_M4, measurevar="Chl.prop", groupvars=c("Site.Geno", "Site", "Genotype"), na.rm=TRUE)

##Plot Average Retention across Treatments
Tol_Chl_M4_SG.plot<-ggplot(Tol_Chl_M4_SG, aes(x=Site.Geno, y=Chl.prop, colour=Genotype)) + 
  geom_errorbar(aes(ymin=Chl.prop-se, ymax=Chl.prop+se), width=cap.sz, linewidth=bar.sz)+
  geom_point(size=point.sz)+ 
   ggtitle("Chlorophyll Retention")+
  theme_classic()+
  theme( axis.title.x = element_text(size = axis.title.sz), 
         axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"),
        axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz),
        legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"),
        legend.position="bottom", 
        legend.direction="horizontal",
        plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
  labs(x="", y="Proportion Retained")+
  ylim(0, 1)+ 
  scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
  scale_color_manual(values = Geno.colors.o)+ 
  stat_pvalue_manual(data=Tol_Chl_M4_lm.p,  y.position=0.55, step.increase=0.20, label="Sig", hide.ns=TRUE); Tol_Chl_M4_SG.plot

Thermal Tolerance Symbionts

Full Set

Run Model

##Check normality
hist(TolData$Sym.prop)

shapiro.test(TolData$Sym.prop)

    Shapiro-Wilk normality test

data:  TolData$Sym.prop
W = 0.92235, p-value = 0.0003707

#Not Normal

##Try log+1 transformation
hist(log(TolData$Sym.prop+1))

shapiro.test(log(TolData$Sym.prop+1))

    Shapiro-Wilk normality test

data:  log(TolData$Sym.prop + 1)
W = 0.92164, p-value = 0.000345

#Not normal 

##Try square transformation
hist((TolData$Sym.prop)^2)

shapiro.test((TolData$Sym.prop)^2)

    Shapiro-Wilk normality test

data:  (TolData$Sym.prop)^2
W = 0.87147, p-value = 3.933e-06

#Not normal 

##Try cubed transformation
hist((TolData$Sym.prop)^3)

shapiro.test((TolData$Sym.prop)^3)

    Shapiro-Wilk normality test

data:  (TolData$Sym.prop)^3
W = 0.80697, p-value = 4.342e-08

#Not normal 


##Model as a function of Site and Genotype and Timepoint
##Model with no transformation and check residuals
Tol_Sym_lm<-lm(Sym.prop~Site+Genotype+TimeP+ Site:Genotype + Site:TimeP + Genotype:TimeP, data=TolData)

Check Residuals

##Check Normality of Residuals
#Distribution 
plot(density(resid(Tol_Sym_lm)))

#Q-Q plot
qqnorm(resid(Tol_Sym_lm)); qqline(resid(Tol_Sym_lm))

##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_Sym_lm), resid(Tol_Sym_lm))

Model Results

#Model Results
summary(Tol_Sym_lm)

Call:
lm(formula = Sym.prop ~ Site + Genotype + TimeP + Site:Genotype + 
    Site:TimeP + Genotype:TimeP, data = TolData)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.37464 -0.08630 -0.00192  0.08972  0.35883 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)         0.544513   0.018576  29.313  < 2e-16 ***
Site.L             -0.001834   0.026197  -0.070 0.944443    
Genotype.L         -0.373179   0.031737 -11.758  < 2e-16 ***
Genotype.Q          0.028101   0.032605   0.862 0.392497    
TimeP.L             0.133503   0.031737   4.206 9.64e-05 ***
TimeP.Q             0.190344   0.032605   5.838 2.94e-07 ***
Site.L:Genotype.L  -0.060389   0.044883  -1.345 0.183999    
Site.L:Genotype.Q   0.074992   0.045930   1.633 0.108240    
Site.L:TimeP.L     -0.159707   0.044883  -3.558 0.000778 ***
Site.L:TimeP.Q     -0.027928   0.045930  -0.608 0.545653    
Genotype.L:TimeP.L  0.085371   0.055316   1.543 0.128489    
Genotype.Q:TimeP.L -0.084924   0.054623  -1.555 0.125749    
Genotype.L:TimeP.Q  0.065986   0.054623   1.208 0.232204    
Genotype.Q:TimeP.Q -0.205029   0.058264  -3.519 0.000878 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1535 on 55 degrees of freedom
Multiple R-squared:  0.8036,    Adjusted R-squared:  0.7572 
F-statistic: 17.31 on 13 and 55 DF,  p-value: 6.022e-15

anova(Tol_Sym_lm)

Analysis of Variance Table

Response: Sym.prop
               Df Sum Sq Mean Sq F value    Pr(>F)    
Site            1 0.0007 0.00068  0.0289  0.865708    
Genotype        2 3.2796 1.63981 69.5880 8.602e-16 ***
TimeP           2 1.1556 0.57780 24.5197 2.438e-08 ***
Site:Genotype   2 0.1111 0.05557  2.3581  0.104096    
Site:TimeP      2 0.3138 0.15688  6.6576  0.002575 ** 
Genotype:TimeP  4 0.4428 0.11069  4.6974  0.002474 ** 
Residuals      55 1.2961 0.02356                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#Effect Size of Predictors
eta_squared(Tol_Sym_lm, partial=FALSE)

# Effect Size for ANOVA (Type I)

Parameter      |     Eta2 |       95% CI
----------------------------------------
Site           | 1.03e-04 | [0.00, 1.00]
Genotype       |     0.50 | [0.33, 1.00]
TimeP          |     0.18 | [0.04, 1.00]
Site:Genotype  |     0.02 | [0.00, 1.00]
Site:TimeP     |     0.05 | [0.00, 1.00]
Genotype:TimeP |     0.07 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
Tol_Sym_lm.res<-data.frame(anova(Tol_Sym_lm))
Tol_Sym_lm.res$Predictor<-rownames(Tol_Sym_lm.res)
Tol_Sym_lm.res$EtaSq<-c(eta_squared(Tol_Sym_lm, partial=FALSE)$Eta2, NA)
Tol_Sym_lm.res$Response<-rep("Symbionts", nrow(Tol_Sym_lm.res))
Tol_Sym_lm.res<-Tol_Sym_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")

Strong influence of Timepoint, including interactions with main variables of interest, Site and Genotype. Will analyze each Timepoint individually.

Sym Summer 1 Timepoint W2

Run Model

##Check normality
hist(TolData_W2$Sym.prop)

shapiro.test(TolData_W2$Sym.prop)

    Shapiro-Wilk normality test

data:  TolData_W2$Sym.prop
W = 0.93605, p-value = 0.1477

#Normal

##Model as a function of Site and Genotype
Tol_Sym_W2_lm<-lm(Sym.prop~Site+Genotype+Site:Genotype, data=TolData_W2)

Check Residuals

##Check Normality of Residuals
#Distribution 
plot(density(resid(Tol_Sym_W2_lm)))

#Q-Q plot
qqnorm(resid(Tol_Sym_W2_lm)); qqline(resid(Tol_Sym_W2_lm))

##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_Sym_W2_lm), resid(Tol_Sym_W2_lm))

Model Results

Overall

#Model Results
summary(Tol_Sym_W2_lm)

Call:
lm(formula = Sym.prop ~ Site + Genotype + Site:Genotype, data = TolData_W2)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.26060 -0.11145  0.00000  0.08525  0.28580 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)        0.5237625  0.0325072  16.112 9.90e-12 ***
Site.L             0.0939568  0.0459722   2.044   0.0568 .  
Genotype.L        -0.4152131  0.0570402  -7.279 1.29e-06 ***
Genotype.Q        -0.0005205  0.0555584  -0.009   0.9926    
Site.L:Genotype.L -0.2061000  0.0806671  -2.555   0.0205 *  
Site.L:Genotype.Q  0.0462891  0.0785715   0.589   0.5635    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.155 on 17 degrees of freedom
Multiple R-squared:  0.7924,    Adjusted R-squared:  0.7313 
F-statistic: 12.98 on 5 and 17 DF,  p-value: 2.647e-05

anova(Tol_Sym_W2_lm)

Analysis of Variance Table

Response: Sym.prop
              Df  Sum Sq Mean Sq F value    Pr(>F)    
Site           1 0.16148 0.16148  6.7208   0.01897 *  
Genotype       2 1.22843 0.61421 25.5640 7.508e-06 ***
Site:Genotype  2 0.16883 0.08441  3.5133   0.05283 .  
Residuals     17 0.40845 0.02403                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#Effect Size of Predictors
eta_squared(Tol_Sym_W2_lm, partial=FALSE)

# Effect Size for ANOVA (Type I)

Parameter     | Eta2 |       95% CI
-----------------------------------
Site          | 0.08 | [0.00, 1.00]
Genotype      | 0.62 | [0.33, 1.00]
Site:Genotype | 0.09 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
Tol_Sym_W2_lm.res<-data.frame(anova(Tol_Sym_W2_lm))
Tol_Sym_W2_lm.res$Predictor<-rownames(Tol_Sym_W2_lm.res)
Tol_Sym_W2_lm.res$EtaSq<-c(eta_squared(Tol_Sym_W2_lm, partial=FALSE)$Eta2, NA)
Tol_Sym_W2_lm.res$Response<-rep("Symbionts", nrow(Tol_Sym_W2_lm.res))
Tol_Sym_W2_lm.res$TimeP<-rep("W2", nrow(Tol_Sym_W2_lm.res))
Tol_Sym_W2_lm.res<-Tol_Sym_W2_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")

Significant effect of Site and Genotype. Marginal effect (p<0.1) of Site * Genotype. Still checking Site*Genotype for comparability across Timepoints.

Pairwise

#Pairwise comparisons across:

#Genotypes within Sites
emmeans(Tol_Sym_W2_lm, pairwise~Genotype | Site)

$emmeans
Site = KL:
 Genotype emmean     SE df lower.CL upper.CL
 AC10      0.634 0.0775 17   0.4708    0.798
 AC12      0.484 0.0775 17   0.3210    0.648
 AC8       0.253 0.0775 17   0.0897    0.417

Site = SS:
 Genotype emmean     SE df lower.CL upper.CL
 AC10      1.000 0.0775 17   0.8365    1.164
 AC12      0.564 0.0775 17   0.4004    0.727
 AC8       0.207 0.0895 17   0.0179    0.396

Confidence level used: 0.95 

$contrasts
Site = KL:
 contrast    estimate    SE df t.ratio p.value
 AC10 - AC12    0.150 0.110 17   1.367  0.3795
 AC10 - AC8     0.381 0.110 17   3.477  0.0077
 AC12 - AC8     0.231 0.110 17   2.110  0.1174

Site = SS:
 contrast    estimate    SE df t.ratio p.value
 AC10 - AC12    0.436 0.110 17   3.979  0.0026
 AC10 - AC8     0.793 0.118 17   6.701  <.0001
 AC12 - AC8     0.357 0.118 17   3.017  0.0201

P value adjustment: tukey method for comparing a family of 3 estimates 

#Sites within Genotypes
emmeans(Tol_Sym_W2_lm, pairwise~Site | Genotype)

$emmeans
Genotype = AC10:
 Site emmean     SE df lower.CL upper.CL
 KL    0.634 0.0775 17   0.4708    0.798
 SS    1.000 0.0775 17   0.8365    1.164

Genotype = AC12:
 Site emmean     SE df lower.CL upper.CL
 KL    0.484 0.0775 17   0.3210    0.648
 SS    0.564 0.0775 17   0.4004    0.727

Genotype = AC8:
 Site emmean     SE df lower.CL upper.CL
 KL    0.253 0.0775 17   0.0897    0.417
 SS    0.207 0.0895 17   0.0179    0.396

Confidence level used: 0.95 

$contrasts
Genotype = AC10:
 contrast estimate    SE df t.ratio p.value
 KL - SS   -0.3657 0.110 17  -3.337  0.0039

Genotype = AC12:
 contrast estimate    SE df t.ratio p.value
 KL - SS   -0.0794 0.110 17  -0.725  0.4785

Genotype = AC8:
 contrast estimate    SE df t.ratio p.value
 KL - SS    0.0465 0.118 17   0.393  0.6994

##Save p-values

#Genotypes within Sites
Tol_Sym_W2_lm.geno<-data.frame(emmeans(Tol_Sym_W2_lm, pairwise~Genotype | Site)$contrasts)
Tol_Sym_W2_lm.geno<-Tol_Sym_W2_lm.geno %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_Sym_W2_lm.geno$group1<-paste(Tol_Sym_W2_lm.geno$Site, Tol_Sym_W2_lm.geno$group1, sep="_")
Tol_Sym_W2_lm.geno$group2<-paste(Tol_Sym_W2_lm.geno$Site, Tol_Sym_W2_lm.geno$group2, sep="_")

#Sites within Genotypes
Tol_Sym_W2_lm.site<-data.frame(emmeans(Tol_Sym_W2_lm, pairwise~Site | Genotype)$contrasts)
Tol_Sym_W2_lm.site<-Tol_Sym_W2_lm.site %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_Sym_W2_lm.site$group1<-paste(Tol_Sym_W2_lm.site$group1, Tol_Sym_W2_lm.site$Genotype, sep="_")
Tol_Sym_W2_lm.site$group2<-paste(Tol_Sym_W2_lm.site$group2, Tol_Sym_W2_lm.site$Genotype, sep="_")

#Full list of p-values
Tol_Sym_W2_lm.p<-rbind(Tol_Sym_W2_lm.geno[,c(1:2,4:8)], Tol_Sym_W2_lm.site[,c(1:2,4:8)])
Tol_Sym_W2_lm.p<-Tol_Sym_W2_lm.p %>% dplyr::rename( p = p.value)

#Add Significance Levels
Tol_Sym_W2_lm.p$Sig<-ifelse(Tol_Sym_W2_lm.p$p<0.001, "***", ifelse(Tol_Sym_W2_lm.p$p<0.01, "**", ifelse(Tol_Sym_W2_lm.p$p<0.05, "*", NA)))

#Specify Response and Timepoint
Tol_Sym_W2_lm.p$Response<-rep("Symbionts", nrow(Tol_Sym_W2_lm.p))
Tol_Sym_W2_lm.p$TimeP<-rep("W2", nrow(Tol_Sym_W2_lm.p))

Plot Retention by Site and Genotype

##Summary statistics by Site and Genotype
Tol_Sym_W2_SG<-summarySE(TolData_W2, measurevar="Sym.prop", groupvars=c("Site.Geno", "Site", "Genotype"), na.rm=TRUE)

##Plot Average Retention across Treatments
Tol_Sym_W2_SG.plot<-ggplot(Tol_Sym_W2_SG, aes(x=Site.Geno, y=Sym.prop, colour=Genotype)) + 
  geom_errorbar(aes(ymin=Sym.prop-se, ymax=Sym.prop+se), width=cap.sz, linewidth=bar.sz)+
  geom_point(size=point.sz)+ 
   ggtitle("Symbiont Retention")+
  theme_classic()+
  theme( axis.title.x = element_text(size = axis.title.sz), 
         axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"),
        axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz),
        legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"),
        legend.position="bottom", 
        legend.direction="horizontal",
        plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
  labs(x="", y="Proportion Retained")+
  ylim(0, 1.5)+ 
  scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
  scale_color_manual(values = Geno.colors.o)+ 
  stat_pvalue_manual(data=Tol_Sym_W2_lm.p,  y.position=0.95, step.increase=0.15, label="Sig", hide.ns=TRUE); Tol_Sym_W2_SG.plot

Sym Summer 2 Timepoint M1

Run Model

##Check normality
hist(TolData_M1$Sym.prop)

shapiro.test(TolData_M1$Sym.prop)

    Shapiro-Wilk normality test

data:  TolData_M1$Sym.prop
W = 0.85748, p-value = 0.004618

#Not normal

##Try log+1 transformation
hist(log(TolData_M1$Sym.prop+1))

shapiro.test(log(TolData_M1$Sym.prop+1))

    Shapiro-Wilk normality test

data:  log(TolData_M1$Sym.prop + 1)
W = 0.87335, p-value = 0.009041

#Not normal but improved

##Model as a function of Site and Genotype
##Model with log transformation
Tol_Sym_M1_lm<-lm(log(Sym.prop+1)~Site+Genotype+Site:Genotype, data=TolData_M1)

Check Residuals

##Check Normality of Residuals
#Distribution 
plot(density(resid(Tol_Sym_M1_lm)))

#Q-Q plot
qqnorm(resid(Tol_Sym_M1_lm)); qqline(resid(Tol_Sym_M1_lm))

##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_Sym_M1_lm), resid(Tol_Sym_M1_lm))

Model Results

Overall

#Model Results
summary(Tol_Sym_M1_lm)

Call:
lm(formula = log(Sym.prop + 1) ~ Site + Genotype + Site:Genotype, 
    data = TolData_M1)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.098079 -0.055799 -0.004888  0.043331  0.151113 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)        0.307716   0.016556  18.587 2.95e-12 ***
Site.L             0.027767   0.023413   1.186 0.252954    
Genotype.L        -0.295494   0.027204 -10.862 8.59e-09 ***
Genotype.Q         0.129387   0.030075   4.302 0.000548 ***
Site.L:Genotype.L  0.080682   0.038472   2.097 0.052224 .  
Site.L:Genotype.Q -0.005572   0.042532  -0.131 0.897412    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.07694 on 16 degrees of freedom
Multiple R-squared:  0.8989,    Adjusted R-squared:  0.8673 
F-statistic: 28.46 on 5 and 16 DF,  p-value: 2e-07

anova(Tol_Sym_M1_lm)

Analysis of Variance Table

Response: log(Sym.prop + 1)
              Df  Sum Sq Mean Sq F value    Pr(>F)    
Site           1 0.00823 0.00823  1.3902    0.2556    
Genotype       2 0.80811 0.40406 68.2486 1.468e-08 ***
Site:Genotype  2 0.02614 0.01307  2.2076    0.1423    
Residuals     16 0.09473 0.00592                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#Effect Size of Predictors
eta_squared(Tol_Sym_M1_lm, partial=FALSE)

# Effect Size for ANOVA (Type I)

Parameter     |     Eta2 |       95% CI
---------------------------------------
Site          | 8.78e-03 | [0.00, 1.00]
Genotype      |     0.86 | [0.73, 1.00]
Site:Genotype |     0.03 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
Tol_Sym_M1_lm.res<-data.frame(anova(Tol_Sym_M1_lm))
Tol_Sym_M1_lm.res$Predictor<-rownames(Tol_Sym_M1_lm.res)
Tol_Sym_M1_lm.res$EtaSq<-c(eta_squared(Tol_Sym_M1_lm, partial=FALSE)$Eta2, NA)
Tol_Sym_M1_lm.res$Response<-rep("Symbionts", nrow(Tol_Sym_M1_lm.res))
Tol_Sym_M1_lm.res$TimeP<-rep("M1", nrow(Tol_Sym_M1_lm.res))
Tol_Sym_M1_lm.res<-Tol_Sym_M1_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")

Significant effect of Genotype. Still checking Site*Genotype for comparability across Timepoints.

Pairwise

#Pairwise comparisons across:

#Genotypes within Sites
emmeans(Tol_Sym_M1_lm, pairwise~Genotype | Site)

$emmeans
Site = KL:
 Genotype emmean     SE df lower.CL upper.CL
 AC10     0.5918 0.0385 16   0.5102    0.673
 AC12     0.1792 0.0444 16   0.0850    0.273
 AC8      0.0932 0.0385 16   0.0117    0.175

Site = SS:
 Genotype emmean     SE df lower.CL upper.CL
 AC10     0.5472 0.0385 16   0.4656    0.629
 AC12     0.2249 0.0444 16   0.1307    0.319
 AC8      0.2100 0.0385 16   0.1284    0.292

Results are given on the log(mu + 1) (not the response) scale. 
Confidence level used: 0.95 

$contrasts
Site = KL:
 contrast    estimate     SE df t.ratio p.value
 AC10 - AC12    0.413 0.0588 16   7.021  <.0001
 AC10 - AC8     0.499 0.0544 16   9.164  <.0001
 AC12 - AC8     0.086 0.0588 16   1.463  0.3339

Site = SS:
 contrast    estimate     SE df t.ratio p.value
 AC10 - AC12    0.322 0.0588 16   5.483  0.0001
 AC10 - AC8     0.337 0.0544 16   6.198  <.0001
 AC12 - AC8     0.015 0.0588 16   0.255  0.9650

Note: contrasts are still on the log(mu + 1) scale 
P value adjustment: tukey method for comparing a family of 3 estimates 

#Sites within Genotypes
emmeans(Tol_Sym_M1_lm, pairwise~Site | Genotype)

$emmeans
Genotype = AC10:
 Site emmean     SE df lower.CL upper.CL
 KL   0.5918 0.0385 16   0.5102    0.673
 SS   0.5472 0.0385 16   0.4656    0.629

Genotype = AC12:
 Site emmean     SE df lower.CL upper.CL
 KL   0.1792 0.0444 16   0.0850    0.273
 SS   0.2249 0.0444 16   0.1307    0.319

Genotype = AC8:
 Site emmean     SE df lower.CL upper.CL
 KL   0.0932 0.0385 16   0.0117    0.175
 SS   0.2100 0.0385 16   0.1284    0.292

Results are given on the log(mu + 1) (not the response) scale. 
Confidence level used: 0.95 

$contrasts
Genotype = AC10:
 contrast estimate     SE df t.ratio p.value
 KL - SS    0.0446 0.0544 16   0.820  0.4241

Genotype = AC12:
 contrast estimate     SE df t.ratio p.value
 KL - SS   -0.0457 0.0628 16  -0.727  0.4775

Genotype = AC8:
 contrast estimate     SE df t.ratio p.value
 KL - SS   -0.1167 0.0544 16  -2.146  0.0476

Note: contrasts are still on the log(mu + 1) scale 

##Save p-values

#Genotypes within Sites
Tol_Sym_M1_lm.geno<-data.frame(emmeans(Tol_Sym_M1_lm, pairwise~Genotype | Site)$contrasts)
Tol_Sym_M1_lm.geno<-Tol_Sym_M1_lm.geno %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_Sym_M1_lm.geno$group1<-paste(Tol_Sym_M1_lm.geno$Site, Tol_Sym_M1_lm.geno$group1, sep="_")
Tol_Sym_M1_lm.geno$group2<-paste(Tol_Sym_M1_lm.geno$Site, Tol_Sym_M1_lm.geno$group2, sep="_")

#Sites within Genotypes
Tol_Sym_M1_lm.site<-data.frame(emmeans(Tol_Sym_M1_lm, pairwise~Site | Genotype)$contrasts)
Tol_Sym_M1_lm.site<-Tol_Sym_M1_lm.site %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_Sym_M1_lm.site$group1<-paste(Tol_Sym_M1_lm.site$group1, Tol_Sym_M1_lm.site$Genotype, sep="_")
Tol_Sym_M1_lm.site$group2<-paste(Tol_Sym_M1_lm.site$group2, Tol_Sym_M1_lm.site$Genotype, sep="_")

#Full list of p-values
Tol_Sym_M1_lm.p<-rbind(Tol_Sym_M1_lm.geno[,c(1:2,4:8)], Tol_Sym_M1_lm.site[,c(1:2,4:8)])
Tol_Sym_M1_lm.p<-Tol_Sym_M1_lm.p %>% dplyr::rename( p = p.value)

#Add Significance Levels
Tol_Sym_M1_lm.p$Sig<-ifelse(Tol_Sym_M1_lm.p$p<0.001, "***", ifelse(Tol_Sym_M1_lm.p$p<0.01, "**", ifelse(Tol_Sym_M1_lm.p$p<0.05, "*", NA)))

#Specify Response and Timepoint
Tol_Sym_M1_lm.p$Response<-rep("Symbionts", nrow(Tol_Sym_M1_lm.p))
Tol_Sym_M1_lm.p$TimeP<-rep("M1", nrow(Tol_Sym_M1_lm.p))

Plot Retention by Site and Genotype

##Summary statistics by Site and Genotype
Tol_Sym_M1_SG<-summarySE(TolData_M1, measurevar="Sym.prop", groupvars=c("Site.Geno", "Site", "Genotype"), na.rm=TRUE)

##Plot Average Retention across Treatments
Tol_Sym_M1_SG.plot<-ggplot(Tol_Sym_M1_SG, aes(x=Site.Geno, y=Sym.prop, colour=Genotype)) + 
  geom_errorbar(aes(ymin=Sym.prop-se, ymax=Sym.prop+se), width=cap.sz, linewidth=bar.sz)+
  geom_point(size=point.sz)+ 
   ggtitle("Symbiont Retention")+
  theme_classic()+
  theme( axis.title.x = element_text(size = axis.title.sz), 
         axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"),
        axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz),
        legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"),
        legend.position="bottom", 
        legend.direction="horizontal",
        plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
  labs(x="", y="Proportion Retained")+
 ylim(0, 1.5)+ 
  scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
  scale_color_manual(values = Geno.colors.o)+ 
  stat_pvalue_manual(data=Tol_Sym_M1_lm.p,  y.position=0.95, step.increase=0.15, label="Sig", hide.ns=TRUE); Tol_Sym_M1_SG.plot

Sym Winter Timepoint M4

Run Model

##Check normality
hist(TolData_M4$Sym.prop)

shapiro.test(TolData_M4$Sym.prop)

    Shapiro-Wilk normality test

data:  TolData_M4$Sym.prop
W = 0.88557, p-value = 0.01078

#Not Normal

##Try square transformation
hist((TolData_M4$Sym.prop)^2)

shapiro.test((TolData_M4$Sym.prop)^2)

    Shapiro-Wilk normality test

data:  (TolData_M4$Sym.prop)^2
W = 0.87973, p-value = 0.008208

#Not Normal

##Try cubed transformation
hist((TolData_M4$Sym.prop)^3)

shapiro.test((TolData_M4$Sym.prop)^3)

    Shapiro-Wilk normality test

data:  (TolData_M4$Sym.prop)^3
W = 0.8488, p-value = 0.002077

#Not Normal

##Model as a function of Site and Genotype
##Model with no transformation and check residuals
Tol_Sym_M4_lm<-lm(Sym.prop~Site+Genotype+Site:Genotype, data=TolData_M4)

Check Residuals

##Check Normality of Residuals
#Distribution 
plot(density(resid(Tol_Sym_M4_lm)))

#Q-Q plot
qqnorm(resid(Tol_Sym_M4_lm)); qqline(resid(Tol_Sym_M4_lm))

##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_Sym_M4_lm), resid(Tol_Sym_M4_lm))

Model Results

Overall

#Model Results
summary(Tol_Sym_M4_lm)

Call:
lm(formula = Sym.prop ~ Site + Genotype + Site:Genotype, data = TolData_M4)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.31820 -0.08804  0.00000  0.08618  0.29740 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)        0.71662    0.03157  22.702 1.07e-14 ***
Site.L            -0.12617    0.04464  -2.826   0.0112 *  
Genotype.L        -0.28587    0.05467  -5.229 5.67e-05 ***
Genotype.Q        -0.11565    0.05467  -2.115   0.0486 *  
Site.L:Genotype.L -0.09494    0.07732  -1.228   0.2353    
Site.L:Genotype.Q  0.17275    0.07732   2.234   0.0384 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1546 on 18 degrees of freedom
Multiple R-squared:  0.7201,    Adjusted R-squared:  0.6423 
F-statistic:  9.26 on 5 and 18 DF,  p-value: 0.0001686

anova(Tol_Sym_M4_lm)

Analysis of Variance Table

Response: Sym.prop
              Df  Sum Sq Mean Sq F value   Pr(>F)    
Site           1 0.19101 0.19101  7.9874 0.011190 *  
Genotype       2 0.76080 0.38040 15.9068 0.000105 ***
Site:Genotype  2 0.15542 0.07771  3.2496 0.062386 .  
Residuals     18 0.43046 0.02391                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#Effect Size of Predictors
eta_squared(Tol_Sym_M4_lm, partial=FALSE)

# Effect Size for ANOVA (Type I)

Parameter     | Eta2 |       95% CI
-----------------------------------
Site          | 0.12 | [0.00, 1.00]
Genotype      | 0.49 | [0.18, 1.00]
Site:Genotype | 0.10 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
Tol_Sym_M4_lm.res<-data.frame(anova(Tol_Sym_M4_lm))
Tol_Sym_M4_lm.res$Predictor<-rownames(Tol_Sym_M4_lm.res)
Tol_Sym_M4_lm.res$EtaSq<-c(eta_squared(Tol_Sym_M4_lm, partial=FALSE)$Eta2, NA)
Tol_Sym_M4_lm.res$Response<-rep("Symbionts", nrow(Tol_Sym_M4_lm.res))
Tol_Sym_M4_lm.res$TimeP<-rep("M4", nrow(Tol_Sym_M4_lm.res))
Tol_Sym_M4_lm.res<-Tol_Sym_M4_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")

Significant effect of Site and Genotype. Marginal effect (p<0.1) of Site * Genotype. Still checking Site*Genotype for comparability across Timepoints.

Pairwise

#Pairwise comparisons across:

#Genotypes within Sites
emmeans(Tol_Sym_M4_lm, pairwise~Genotype | Site)

$emmeans
Site = KL:
 Genotype emmean     SE df lower.CL upper.CL
 AC10      0.863 0.0773 18    0.701    1.026
 AC12      1.000 0.0773 18    0.838    1.162
 AC8       0.554 0.0773 18    0.392    0.717

Site = SS:
 Genotype emmean     SE df lower.CL upper.CL
 AC10      0.880 0.0773 18    0.717    1.042
 AC12      0.622 0.0773 18    0.460    0.785
 AC8       0.380 0.0773 18    0.218    0.543

Confidence level used: 0.95 

$contrasts
Site = KL:
 contrast    estimate    SE df t.ratio p.value
 AC10 - AC12   -0.137 0.109 18  -1.249  0.4409
 AC10 - AC8     0.309 0.109 18   2.829  0.0285
 AC12 - AC8     0.446 0.109 18   4.078  0.0019

Site = SS:
 contrast    estimate    SE df t.ratio p.value
 AC10 - AC12    0.258 0.109 18   2.356  0.0734
 AC10 - AC8     0.499 0.109 18   4.565  0.0007
 AC12 - AC8     0.242 0.109 18   2.210  0.0965

P value adjustment: tukey method for comparing a family of 3 estimates 

#Sites within Genotypes
emmeans(Tol_Sym_M4_lm, pairwise~Site | Genotype)

$emmeans
Genotype = AC10:
 Site emmean     SE df lower.CL upper.CL
 KL    0.863 0.0773 18    0.701    1.026
 SS    0.880 0.0773 18    0.717    1.042

Genotype = AC12:
 Site emmean     SE df lower.CL upper.CL
 KL    1.000 0.0773 18    0.838    1.162
 SS    0.622 0.0773 18    0.460    0.785

Genotype = AC8:
 Site emmean     SE df lower.CL upper.CL
 KL    0.554 0.0773 18    0.392    0.717
 SS    0.380 0.0773 18    0.218    0.543

Confidence level used: 0.95 

$contrasts
Genotype = AC10:
 contrast estimate    SE df t.ratio p.value
 KL - SS   -0.0163 0.109 18  -0.149  0.8835

Genotype = AC12:
 contrast estimate    SE df t.ratio p.value
 KL - SS    0.3779 0.109 18   3.456  0.0028

Genotype = AC8:
 contrast estimate    SE df t.ratio p.value
 KL - SS    0.1736 0.109 18   1.588  0.1297

##Save p-values

#Genotypes within Sites
Tol_Sym_M4_lm.geno<-data.frame(emmeans(Tol_Sym_M4_lm, pairwise~Genotype | Site)$contrasts)
Tol_Sym_M4_lm.geno<-Tol_Sym_M4_lm.geno %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_Sym_M4_lm.geno$group1<-paste(Tol_Sym_M4_lm.geno$Site, Tol_Sym_M4_lm.geno$group1, sep="_")
Tol_Sym_M4_lm.geno$group2<-paste(Tol_Sym_M4_lm.geno$Site, Tol_Sym_M4_lm.geno$group2, sep="_")

#Sites within Genotypes
Tol_Sym_M4_lm.site<-data.frame(emmeans(Tol_Sym_M4_lm, pairwise~Site | Genotype)$contrasts)
Tol_Sym_M4_lm.site<-Tol_Sym_M4_lm.site %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_Sym_M4_lm.site$group1<-paste(Tol_Sym_M4_lm.site$group1, Tol_Sym_M4_lm.site$Genotype, sep="_")
Tol_Sym_M4_lm.site$group2<-paste(Tol_Sym_M4_lm.site$group2, Tol_Sym_M4_lm.site$Genotype, sep="_")

#Full list of p-values
Tol_Sym_M4_lm.p<-rbind(Tol_Sym_M4_lm.geno[,c(1:2,4:8)], Tol_Sym_M4_lm.site[,c(1:2,4:8)])
Tol_Sym_M4_lm.p<-Tol_Sym_M4_lm.p %>% dplyr::rename( p = p.value)

#Add Significance Levels
Tol_Sym_M4_lm.p$Sig<-ifelse(Tol_Sym_M4_lm.p$p<0.001, "***", ifelse(Tol_Sym_M4_lm.p$p<0.01, "**", ifelse(Tol_Sym_M4_lm.p$p<0.05, "*", NA)))

#Specify Response and Timepoint
Tol_Sym_M4_lm.p$Response<-rep("Symbionts", nrow(Tol_Sym_M4_lm.p))
Tol_Sym_M4_lm.p$TimeP<-rep("M4", nrow(Tol_Sym_M4_lm.p))

Plot Retention by Site and Genotype

##Summary statistics by Site and Genotype
Tol_Sym_M4_SG<-summarySE(TolData_M4, measurevar="Sym.prop", groupvars=c("Site.Geno", "Site", "Genotype"), na.rm=TRUE)

##Plot Average Retention across Treatments
Tol_Sym_M4_SG.plot<-ggplot(Tol_Sym_M4_SG, aes(x=Site.Geno, y=Sym.prop, colour=Genotype)) + 
  geom_errorbar(aes(ymin=Sym.prop-se, ymax=Sym.prop+se), width=cap.sz, linewidth=bar.sz)+
  geom_point(size=point.sz)+ 
   ggtitle("Symbiont Retention")+
  theme_classic()+
  theme( axis.title.x = element_text(size = axis.title.sz), 
         axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"),
        axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz),
        legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"),
        legend.position="bottom", 
        legend.direction="horizontal",
        plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
  labs(x="", y="Proportion Retained")+
  ylim(0, 1.5)+ 
  scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
  scale_color_manual(values = Geno.colors.o)+ 
  stat_pvalue_manual(data=Tol_Sym_M4_lm.p,  y.position=1.1, step.increase=0.15, label="Sig", hide.ns=TRUE); Tol_Sym_M4_SG.plot

Thermal Tolerance Color Full Set

Full Set

Run Model

##Check normality
hist(TolData$Score_Full.prop)

shapiro.test(TolData$Score_Full.prop)

    Shapiro-Wilk normality test

data:  TolData$Score_Full.prop
W = 0.90488, p-value = 6.821e-05

#Not Normal

##Try square transformation
hist((TolData$Score_Full.prop)^2)

shapiro.test((TolData$Score_Full.prop)^2)

    Shapiro-Wilk normality test

data:  (TolData$Score_Full.prop)^2
W = 0.9194, p-value = 0.0002754

#Not normal 

##Try cubed transformation
hist((TolData$Score_Full.prop)^3)

shapiro.test((TolData$Score_Full.prop)^3)

    Shapiro-Wilk normality test

data:  (TolData$Score_Full.prop)^3
W = 0.91645, p-value = 0.0002055

#Not normal 


##Model as a function of Site and Genotype and Timepoint
##Model with no transformation and check residuals
Tol_ScoreF_lm<-lm(Score_Full.prop~Site+Genotype+TimeP+ Site:Genotype + Site:TimeP + Genotype:TimeP, data=TolData)

Check Residuals

##Check Normality of Residuals
#Distribution 
plot(density(resid(Tol_ScoreF_lm)))

#Q-Q plot
qqnorm(resid(Tol_ScoreF_lm)); qqline(resid(Tol_ScoreF_lm))

##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_ScoreF_lm), resid(Tol_ScoreF_lm))

Residuals are not great, need to check for other modeling options for Color Full Set.

Model Results

#Model Results
summary(Tol_ScoreF_lm)

Call:
lm(formula = Score_Full.prop ~ Site + Genotype + TimeP + Site:Genotype + 
    Site:TimeP + Genotype:TimeP, data = TolData)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.217253 -0.020809  0.001419  0.037990  0.181147 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)         0.816721   0.010646  76.716  < 2e-16 ***
Site.L              0.068538   0.015014   4.565 2.85e-05 ***
Genotype.L         -0.162107   0.018189  -8.912 2.92e-12 ***
Genotype.Q         -0.000730   0.018686  -0.039  0.96898    
TimeP.L             0.143092   0.018189   7.867 1.44e-10 ***
TimeP.Q            -0.021653   0.018686  -1.159  0.25157    
Site.L:Genotype.L   0.033429   0.025724   1.300  0.19918    
Site.L:Genotype.Q  -0.007619   0.026324  -0.289  0.77333    
Site.L:TimeP.L     -0.039262   0.025724  -1.526  0.13267    
Site.L:TimeP.Q     -0.005353   0.026324  -0.203  0.83961    
Genotype.L:TimeP.L  0.055269   0.031703   1.743  0.08686 .  
Genotype.Q:TimeP.L  0.032487   0.031306   1.038  0.30393    
Genotype.L:TimeP.Q -0.003858   0.031306  -0.123  0.90237    
Genotype.Q:TimeP.Q -0.108144   0.033392  -3.239  0.00204 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.08798 on 55 degrees of freedom
Multiple R-squared:  0.7679,    Adjusted R-squared:  0.7131 
F-statistic:    14 on 13 and 55 DF,  p-value: 4.691e-13

anova(Tol_ScoreF_lm)

Analysis of Variance Table

Response: Score_Full.prop
               Df  Sum Sq  Mean Sq F value    Pr(>F)    
Site            1 0.18599 0.185989 24.0291 8.773e-06 ***
Genotype        2 0.58580 0.292900 37.8415 4.613e-11 ***
TimeP           2 0.48992 0.244959 31.6477 7.135e-10 ***
Site:Genotype   2 0.01497 0.007484  0.9669   0.38664    
Site:TimeP      2 0.02075 0.010376  1.3405   0.27013    
Genotype:TimeP  4 0.11112 0.027780  3.5891   0.01137 *  
Residuals      55 0.42571 0.007740                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#Effect Size of Predictors
eta_squared(Tol_ScoreF_lm, partial=FALSE)

# Effect Size for ANOVA (Type I)

Parameter      |     Eta2 |       95% CI
----------------------------------------
Site           |     0.10 | [0.01, 1.00]
Genotype       |     0.32 | [0.15, 1.00]
TimeP          |     0.27 | [0.10, 1.00]
Site:Genotype  | 8.16e-03 | [0.00, 1.00]
Site:TimeP     |     0.01 | [0.00, 1.00]
Genotype:TimeP |     0.06 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
Tol_ScoreF_lm.res<-data.frame(anova(Tol_ScoreF_lm))
Tol_ScoreF_lm.res$Predictor<-rownames(Tol_ScoreF_lm.res)
Tol_ScoreF_lm.res$EtaSq<-c(eta_squared(Tol_ScoreF_lm, partial=FALSE)$Eta2, NA)
Tol_ScoreF_lm.res$Response<-rep("Color_FullSet", nrow(Tol_ScoreF_lm.res))
Tol_ScoreF_lm.res<-Tol_ScoreF_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")

Strong influence of Timepoint, including interactions with Genotype. Will analyze each Timepoint individually.

Color Score Full Set Summer 1 Timepoint W2

Run Model

##Check normality
hist(TolData_W2$Score_Full.prop)

shapiro.test(TolData_W2$Score_Full.prop)

    Shapiro-Wilk normality test

data:  TolData_W2$Score_Full.prop
W = 0.96098, p-value = 0.4834

#Normal

##Model as a function of Site and Genotype
Tol_ScoreF_W2_lm<-lm(Score_Full.prop~Site+Genotype+Site:Genotype, data=TolData_W2)

Check Residuals

##Check Normality of Residuals
#Distribution 
plot(density(resid(Tol_ScoreF_W2_lm)))

#Q-Q plot
qqnorm(resid(Tol_ScoreF_W2_lm)); qqline(resid(Tol_ScoreF_W2_lm))

##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_ScoreF_W2_lm), resid(Tol_ScoreF_W2_lm))

Model Results

Overall

#Model Results
summary(Tol_ScoreF_W2_lm)

Call:
lm(formula = Score_Full.prop ~ Site + Genotype + Site:Genotype, 
    data = TolData_W2)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.20185 -0.04952 -0.00045  0.04135  0.16725 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)        0.706208   0.023330  30.270 3.15e-16 ***
Site.L             0.093421   0.032994   2.831 0.011517 *  
Genotype.L        -0.203806   0.040938  -4.978 0.000115 ***
Genotype.Q        -0.068453   0.039874  -1.717 0.104193    
Site.L:Genotype.L  0.003375   0.057895   0.058 0.954193    
Site.L:Genotype.Q  0.010407   0.056391   0.185 0.855768    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1112 on 17 degrees of freedom
Multiple R-squared:  0.6882,    Adjusted R-squared:  0.5964 
F-statistic: 7.503 on 5 and 17 DF,  p-value: 7e-04

anova(Tol_ScoreF_W2_lm)

Analysis of Variance Table

Response: Score_Full.prop
              Df  Sum Sq  Mean Sq F value    Pr(>F)    
Site           1 0.12458 0.124577 10.0661 0.0055647 ** 
Genotype       2 0.33925 0.169625 13.7060 0.0002851 ***
Site:Genotype  2 0.00045 0.000226  0.0183 0.9819098    
Residuals     17 0.21039 0.012376                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#Effect Size of Predictors
eta_squared(Tol_ScoreF_W2_lm, partial=FALSE)

# Effect Size for ANOVA (Type I)

Parameter     |     Eta2 |       95% CI
---------------------------------------
Site          |     0.18 | [0.00, 1.00]
Genotype      |     0.50 | [0.17, 1.00]
Site:Genotype | 6.70e-04 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
Tol_ScoreF_W2_lm.res<-data.frame(anova(Tol_ScoreF_W2_lm))
Tol_ScoreF_W2_lm.res$Predictor<-rownames(Tol_ScoreF_W2_lm.res)
Tol_ScoreF_W2_lm.res$EtaSq<-c(eta_squared(Tol_ScoreF_W2_lm, partial=FALSE)$Eta2, NA)
Tol_ScoreF_W2_lm.res$Response<-rep("Color_FullSet", nrow(Tol_ScoreF_W2_lm.res))
Tol_ScoreF_W2_lm.res$TimeP<-rep("W2", nrow(Tol_ScoreF_W2_lm.res))
Tol_ScoreF_W2_lm.res<-Tol_ScoreF_W2_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")

Significant effect of Site and Genotype. Still checking Site*Genotype for comparability across Timepoints.

Pairwise

#Pairwise comparisons across:

#Genotypes within Sites
emmeans(Tol_ScoreF_W2_lm, pairwise~Genotype | Site)

$emmeans
Site = KL:
 Genotype emmean     SE df lower.CL upper.CL
 AC10      0.755 0.0556 17    0.638    0.872
 AC12      0.702 0.0556 17    0.585    0.819
 AC8       0.463 0.0556 17    0.346    0.581

Site = SS:
 Genotype emmean     SE df lower.CL upper.CL
 AC10      0.890 0.0556 17    0.772    1.007
 AC12      0.822 0.0556 17    0.705    0.940
 AC8       0.605 0.0642 17    0.469    0.740

Confidence level used: 0.95 

$contrasts
Site = KL:
 contrast    estimate     SE df t.ratio p.value
 AC10 - AC12   0.0529 0.0787 17   0.673  0.7820
 AC10 - AC8    0.2916 0.0787 17   3.707  0.0047
 AC12 - AC8    0.2387 0.0787 17   3.034  0.0195

Site = SS:
 contrast    estimate     SE df t.ratio p.value
 AC10 - AC12   0.0676 0.0787 17   0.859  0.6723
 AC10 - AC8    0.2848 0.0850 17   3.353  0.0100
 AC12 - AC8    0.2172 0.0850 17   2.557  0.0508

P value adjustment: tukey method for comparing a family of 3 estimates 

#Sites within Genotypes
emmeans(Tol_ScoreF_W2_lm, pairwise~Site | Genotype)

$emmeans
Genotype = AC10:
 Site emmean     SE df lower.CL upper.CL
 KL    0.755 0.0556 17    0.638    0.872
 SS    0.890 0.0556 17    0.772    1.007

Genotype = AC12:
 Site emmean     SE df lower.CL upper.CL
 KL    0.702 0.0556 17    0.585    0.819
 SS    0.822 0.0556 17    0.705    0.940

Genotype = AC8:
 Site emmean     SE df lower.CL upper.CL
 KL    0.463 0.0556 17    0.346    0.581
 SS    0.605 0.0642 17    0.469    0.740

Confidence level used: 0.95 

$contrasts
Genotype = AC10:
 contrast estimate     SE df t.ratio p.value
 KL - SS    -0.135 0.0787 17  -1.713  0.1049

Genotype = AC12:
 contrast estimate     SE df t.ratio p.value
 KL - SS    -0.120 0.0787 17  -1.527  0.1452

Genotype = AC8:
 contrast estimate     SE df t.ratio p.value
 KL - SS    -0.141 0.0850 17  -1.665  0.1142

##Save p-values

#Genotypes within Sites
Tol_ScoreF_W2_lm.geno<-data.frame(emmeans(Tol_ScoreF_W2_lm, pairwise~Genotype | Site)$contrasts)
Tol_ScoreF_W2_lm.geno<-Tol_ScoreF_W2_lm.geno %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_ScoreF_W2_lm.geno$group1<-paste(Tol_ScoreF_W2_lm.geno$Site, Tol_ScoreF_W2_lm.geno$group1, sep="_")
Tol_ScoreF_W2_lm.geno$group2<-paste(Tol_ScoreF_W2_lm.geno$Site, Tol_ScoreF_W2_lm.geno$group2, sep="_")

#Sites within Genotypes
Tol_ScoreF_W2_lm.site<-data.frame(emmeans(Tol_ScoreF_W2_lm, pairwise~Site | Genotype)$contrasts)
Tol_ScoreF_W2_lm.site<-Tol_ScoreF_W2_lm.site %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_ScoreF_W2_lm.site$group1<-paste(Tol_ScoreF_W2_lm.site$group1, Tol_ScoreF_W2_lm.site$Genotype, sep="_")
Tol_ScoreF_W2_lm.site$group2<-paste(Tol_ScoreF_W2_lm.site$group2, Tol_ScoreF_W2_lm.site$Genotype, sep="_")

#Full list of p-values
Tol_ScoreF_W2_lm.p<-rbind(Tol_ScoreF_W2_lm.geno[,c(1:2,4:8)], Tol_ScoreF_W2_lm.site[,c(1:2,4:8)])
Tol_ScoreF_W2_lm.p<-Tol_ScoreF_W2_lm.p %>% dplyr::rename( p = p.value)

#Add Significance Levels
Tol_ScoreF_W2_lm.p$Sig<-ifelse(Tol_ScoreF_W2_lm.p$p<0.001, "***", ifelse(Tol_ScoreF_W2_lm.p$p<0.01, "**", ifelse(Tol_ScoreF_W2_lm.p$p<0.05, "*", NA)))

#Specify Response and Timepoint
Tol_ScoreF_W2_lm.p$Response<-rep("Color_FullSet", nrow(Tol_ScoreF_W2_lm.p))
Tol_ScoreF_W2_lm.p$TimeP<-rep("W2", nrow(Tol_ScoreF_W2_lm.p))

Plot Retention by Site and Genotype

##Summary statistics by Site and Genotype
Tol_ScoreF_W2_SG<-summarySE(TolData_W2, measurevar="Score_Full.prop", groupvars=c("Site.Geno", "Site", "Genotype"), na.rm=TRUE)

##Plot Average Retention across Treatments
Tol_ScoreF_W2_SG.plot<-ggplot(Tol_ScoreF_W2_SG, aes(x=Site.Geno, y=Score_Full.prop, colour=Genotype)) + 
  geom_errorbar(aes(ymin=Score_Full.prop-se, ymax=Score_Full.prop+se), width=cap.sz, linewidth=bar.sz)+
  geom_point(size=point.sz)+ 
   ggtitle("Color Retention Full Set")+
  theme_classic()+
  theme( axis.title.x = element_text(size = axis.title.sz), 
         axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"),
        axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz),
        legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"),
        legend.position="bottom", 
        legend.direction="horizontal",
        plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
  labs(x="", y="Proportion Retained")+
  ylim(0, 1.5)+ 
  scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
  scale_color_manual(values = Geno.colors.o)+ 
  stat_pvalue_manual(data=Tol_ScoreF_W2_lm.p,  y.position=0.9, step.increase=0.35, label="Sig", hide.ns=TRUE); Tol_ScoreF_W2_SG.plot

Color Score Full Set Summer 2 Timepoint M1

Run Model

##Check normality
hist(TolData_M1$Score_Full.prop)

shapiro.test(TolData_M1$Score_Full.prop)

    Shapiro-Wilk normality test

data:  TolData_M1$Score_Full.prop
W = 0.89042, p-value = 0.01921

#Not normal

##Try square transformation
hist((TolData_M1$Score_Full.prop)^2)

shapiro.test((TolData_M1$Score_Full.prop)^2)

    Shapiro-Wilk normality test

data:  (TolData_M1$Score_Full.prop)^2
W = 0.89427, p-value = 0.02287

#Not normal

##Try cubed transformation
hist((TolData_M1$Score_Full.prop)^3)

shapiro.test((TolData_M1$Score_Full.prop)^3)

    Shapiro-Wilk normality test

data:  (TolData_M1$Score_Full.prop)^3
W = 0.88949, p-value = 0.01842

#Not normal

##Model as a function of Site and Genotype
##Model with no transformation and check residuals
Tol_ScoreF_M1_lm<-lm(Score_Full.prop~Site+Genotype+Site:Genotype, data=TolData_M1)

Check Residuals

##Check Normality of Residuals
#Distribution 
plot(density(resid(Tol_ScoreF_M1_lm)))

#Q-Q plot
qqnorm(resid(Tol_ScoreF_M1_lm)); qqline(resid(Tol_ScoreF_M1_lm))

##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_ScoreF_M1_lm), resid(Tol_ScoreF_M1_lm))

Model Results

Overall

#Model Results
summary(Tol_ScoreF_M1_lm)

Call:
lm(formula = Score_Full.prop ~ Site + Genotype + Site:Genotype, 
    data = TolData_M1)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.205567 -0.027956  0.001925  0.027094  0.192833 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)        0.83440    0.02043  40.846  < 2e-16 ***
Site.L             0.07153    0.02889   2.476 0.024843 *  
Genotype.L        -0.15896    0.03357  -4.736 0.000224 ***
Genotype.Q         0.08757    0.03711   2.360 0.031325 *  
Site.L:Genotype.L  0.07642    0.04747   1.610 0.126961    
Site.L:Genotype.Q  0.01094    0.05248   0.208 0.837562    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.09494 on 16 degrees of freedom
Multiple R-squared:  0.6982,    Adjusted R-squared:  0.6039 
F-statistic: 7.403 on 5 and 16 DF,  p-value: 0.0009114

anova(Tol_ScoreF_M1_lm)

Analysis of Variance Table

Response: Score_Full.prop
              Df   Sum Sq  Mean Sq F value    Pr(>F)    
Site           1 0.057569 0.057569  6.3867 0.0224118 *  
Genotype       2 0.252333 0.126167 13.9969 0.0003061 ***
Site:Genotype  2 0.023755 0.011877  1.3177 0.2952999    
Residuals     16 0.144222 0.009014                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#Effect Size of Predictors
eta_squared(Tol_ScoreF_M1_lm, partial=FALSE)

# Effect Size for ANOVA (Type I)

Parameter     | Eta2 |       95% CI
-----------------------------------
Site          | 0.12 | [0.00, 1.00]
Genotype      | 0.53 | [0.19, 1.00]
Site:Genotype | 0.05 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
Tol_ScoreF_M1_lm.res<-data.frame(anova(Tol_ScoreF_M1_lm))
Tol_ScoreF_M1_lm.res$Predictor<-rownames(Tol_ScoreF_M1_lm.res)
Tol_ScoreF_M1_lm.res$EtaSq<-c(eta_squared(Tol_ScoreF_M1_lm, partial=FALSE)$Eta2, NA)
Tol_ScoreF_M1_lm.res$Response<-rep("Color_FullSet", nrow(Tol_ScoreF_M1_lm.res))
Tol_ScoreF_M1_lm.res$TimeP<-rep("M1", nrow(Tol_ScoreF_M1_lm.res))
Tol_ScoreF_M1_lm.res<-Tol_ScoreF_M1_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")

Significant effect of Site and Genotype. Still checking Site*Genotype for comparability across Timepoints.

Pairwise

#Pairwise comparisons across:

#Genotypes within Sites
emmeans(Tol_ScoreF_M1_lm, pairwise~Genotype | Site)

$emmeans
Site = KL:
 Genotype emmean     SE df lower.CL upper.CL
 AC10      0.967 0.0475 16    0.866    1.068
 AC12      0.719 0.0548 16    0.602    0.835
 AC8       0.666 0.0475 16    0.565    0.766

Site = SS:
 Genotype emmean     SE df lower.CL upper.CL
 AC10      0.998 0.0475 16    0.897    1.099
 AC12      0.807 0.0548 16    0.691    0.923
 AC8       0.850 0.0475 16    0.749    0.950

Confidence level used: 0.95 

$contrasts
Site = KL:
 contrast    estimate     SE df t.ratio p.value
 AC10 - AC12   0.2484 0.0725 16   3.425  0.0092
 AC10 - AC8    0.3012 0.0671 16   4.487  0.0010
 AC12 - AC8    0.0528 0.0725 16   0.729  0.7504

Site = SS:
 contrast    estimate     SE df t.ratio p.value
 AC10 - AC12   0.1909 0.0725 16   2.633  0.0452
 AC10 - AC8    0.1484 0.0671 16   2.210  0.0998
 AC12 - AC8   -0.0425 0.0725 16  -0.587  0.8292

P value adjustment: tukey method for comparing a family of 3 estimates 

#Sites within Genotypes
emmeans(Tol_ScoreF_M1_lm, pairwise~Site | Genotype)

$emmeans
Genotype = AC10:
 Site emmean     SE df lower.CL upper.CL
 KL    0.967 0.0475 16    0.866    1.068
 SS    0.998 0.0475 16    0.897    1.099

Genotype = AC12:
 Site emmean     SE df lower.CL upper.CL
 KL    0.719 0.0548 16    0.602    0.835
 SS    0.807 0.0548 16    0.691    0.923

Genotype = AC8:
 Site emmean     SE df lower.CL upper.CL
 KL    0.666 0.0475 16    0.565    0.766
 SS    0.850 0.0475 16    0.749    0.950

Confidence level used: 0.95 

$contrasts
Genotype = AC10:
 contrast estimate     SE df t.ratio p.value
 KL - SS   -0.0311 0.0671 16  -0.463  0.6499

Genotype = AC12:
 contrast estimate     SE df t.ratio p.value
 KL - SS   -0.0885 0.0775 16  -1.142  0.2702

Genotype = AC8:
 contrast estimate     SE df t.ratio p.value
 KL - SS   -0.1839 0.0671 16  -2.739  0.0146

##Save p-values

#Genotypes within Sites
Tol_ScoreF_M1_lm.geno<-data.frame(emmeans(Tol_ScoreF_M1_lm, pairwise~Genotype | Site)$contrasts)
Tol_ScoreF_M1_lm.geno<-Tol_ScoreF_M1_lm.geno %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_ScoreF_M1_lm.geno$group1<-paste(Tol_ScoreF_M1_lm.geno$Site, Tol_ScoreF_M1_lm.geno$group1, sep="_")
Tol_ScoreF_M1_lm.geno$group2<-paste(Tol_ScoreF_M1_lm.geno$Site, Tol_ScoreF_M1_lm.geno$group2, sep="_")

#Sites within Genotypes
Tol_ScoreF_M1_lm.site<-data.frame(emmeans(Tol_ScoreF_M1_lm, pairwise~Site | Genotype)$contrasts)
Tol_ScoreF_M1_lm.site<-Tol_ScoreF_M1_lm.site %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_ScoreF_M1_lm.site$group1<-paste(Tol_ScoreF_M1_lm.site$group1, Tol_ScoreF_M1_lm.site$Genotype, sep="_")
Tol_ScoreF_M1_lm.site$group2<-paste(Tol_ScoreF_M1_lm.site$group2, Tol_ScoreF_M1_lm.site$Genotype, sep="_")

#Full list of p-values
Tol_ScoreF_M1_lm.p<-rbind(Tol_ScoreF_M1_lm.geno[,c(1:2,4:8)], Tol_ScoreF_M1_lm.site[,c(1:2,4:8)])
Tol_ScoreF_M1_lm.p<-Tol_ScoreF_M1_lm.p %>% dplyr::rename( p = p.value)

#Add Significance Levels
Tol_ScoreF_M1_lm.p$Sig<-ifelse(Tol_ScoreF_M1_lm.p$p<0.001, "***", ifelse(Tol_ScoreF_M1_lm.p$p<0.01, "**", ifelse(Tol_ScoreF_M1_lm.p$p<0.05, "*", NA)))

#Specify Response and Timepoint
Tol_ScoreF_M1_lm.p$Response<-rep("Color_FullSet", nrow(Tol_ScoreF_M1_lm.p))
Tol_ScoreF_M1_lm.p$TimeP<-rep("M1", nrow(Tol_ScoreF_M1_lm.p))

Plot Retention by Site and Genotype

##Summary statistics by Site and Genotype
Tol_ScoreF_M1_SG<-summarySE(TolData_M1, measurevar="Score_Full.prop", groupvars=c("Site.Geno", "Site", "Genotype"), na.rm=TRUE)

##Plot Average Retention across Treatments
Tol_ScoreF_M1_SG.plot<-ggplot(Tol_ScoreF_M1_SG, aes(x=Site.Geno, y=Score_Full.prop, colour=Genotype)) + 
  geom_errorbar(aes(ymin=Score_Full.prop-se, ymax=Score_Full.prop+se), width=cap.sz, linewidth=bar.sz)+
  geom_point(size=point.sz)+ 
   ggtitle("Color Retention Full Set")+
  theme_classic()+
  theme( axis.title.x = element_text(size = axis.title.sz), 
         axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"),
        axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz),
        legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"),
        legend.position="bottom", 
        legend.direction="horizontal",
        plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
  labs(x="", y="Proportion Retained")+
 ylim(0, 1.5)+ 
  scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
  scale_color_manual(values = Geno.colors.o)+ 
  stat_pvalue_manual(data=Tol_ScoreF_M1_lm.p,  y.position=1.1, step.increase=0.2, label="Sig", hide.ns=TRUE); Tol_ScoreF_M1_SG.plot

Color Score Full Set Winter Timepoint M4

Run Model

##Check normality
hist(TolData_M4$Score_Full.prop)

shapiro.test(TolData_M4$Score_Full.prop)

    Shapiro-Wilk normality test

data:  TolData_M4$Score_Full.prop
W = 0.83109, p-value = 0.0009936

#Not Normal

##Try square transformation
hist((TolData_M4$Score_Full.prop)^2)

shapiro.test((TolData_M4$Score_Full.prop)^2)

    Shapiro-Wilk normality test

data:  (TolData_M4$Score_Full.prop)^2
W = 0.8317, p-value = 0.001019

#Not Normal

##Try cubed transformation
hist((TolData_M4$Score_Full.prop)^3)

shapiro.test((TolData_M4$Score_Full.prop)^3)

    Shapiro-Wilk normality test

data:  (TolData_M4$Score_Full.prop)^3
W = 0.83161, p-value = 0.001015

#Not Normal

##Model as a function of Site and Genotype
##Model with no transformation and check residuals
Tol_ScoreF_M4_lm<-lm(Score_Full.prop~Site+Genotype+Site:Genotype, data=TolData_M4)

Check Residuals

##Check Normality of Residuals
#Distribution 
plot(density(resid(Tol_ScoreF_M4_lm)))

#Q-Q plot
qqnorm(resid(Tol_ScoreF_M4_lm)); qqline(resid(Tol_ScoreF_M4_lm))

##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_ScoreF_M4_lm), resid(Tol_ScoreF_M4_lm))

Model Results

Overall

#Model Results
summary(Tol_ScoreF_M4_lm)

Call:
lm(formula = Score_Full.prop ~ Site + Genotype + Site:Genotype, 
    data = TolData_M4)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.137800 -0.019025  0.003675  0.020575  0.099500 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)        0.90906    0.01100  82.614  < 2e-16 ***
Site.L             0.03859    0.01556   2.480   0.0233 *  
Genotype.L        -0.12460    0.01906  -6.538 3.82e-06 ***
Genotype.Q        -0.02191    0.01906  -1.149   0.2654    
Site.L:Genotype.L  0.01901    0.02695   0.705   0.4896    
Site.L:Genotype.Q -0.04168    0.02695  -1.546   0.1394    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.05391 on 18 degrees of freedom
Multiple R-squared:  0.7468,    Adjusted R-squared:  0.6765 
F-statistic: 10.62 on 5 and 18 DF,  p-value: 7.153e-05

anova(Tol_ScoreF_M4_lm)

Analysis of Variance Table

Response: Score_Full.prop
              Df   Sum Sq  Mean Sq F value    Pr(>F)    
Site           1 0.017871 0.017871  6.1496   0.02326 *  
Genotype       2 0.128043 0.064021 22.0308 1.452e-05 ***
Site:Genotype  2 0.008394 0.004197  1.4442   0.26199    
Residuals     18 0.052308 0.002906                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#Effect Size of Predictors
eta_squared(Tol_ScoreF_M4_lm, partial=FALSE)

# Effect Size for ANOVA (Type I)

Parameter     | Eta2 |       95% CI
-----------------------------------
Site          | 0.09 | [0.00, 1.00]
Genotype      | 0.62 | [0.34, 1.00]
Site:Genotype | 0.04 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
Tol_ScoreF_M4_lm.res<-data.frame(anova(Tol_ScoreF_M4_lm))
Tol_ScoreF_M4_lm.res$Predictor<-rownames(Tol_ScoreF_M4_lm.res)
Tol_ScoreF_M4_lm.res$EtaSq<-c(eta_squared(Tol_ScoreF_M4_lm, partial=FALSE)$Eta2, NA)
Tol_ScoreF_M4_lm.res$Response<-rep("Color_FullSet", nrow(Tol_ScoreF_M4_lm.res))
Tol_ScoreF_M4_lm.res$TimeP<-rep("M4", nrow(Tol_ScoreF_M4_lm.res))
Tol_ScoreF_M4_lm.res<-Tol_ScoreF_M4_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")

Significant effect of Site and Genotype. Still checking Site*Genotype for comparability across Timepoints.

Pairwise

#Pairwise comparisons across:

#Genotypes within Sites
emmeans(Tol_ScoreF_M4_lm, pairwise~Genotype | Site)

$emmeans
Site = KL:
 Genotype emmean    SE df lower.CL upper.CL
 AC10      0.982 0.027 18    0.926    1.039
 AC12      0.876 0.027 18    0.819    0.932
 AC8       0.787 0.027 18    0.731    0.844

Site = SS:
 Genotype emmean    SE df lower.CL upper.CL
 AC10      0.994 0.027 18    0.937    1.051
 AC12      0.978 0.027 18    0.922    1.035
 AC8       0.837 0.027 18    0.780    0.893

Confidence level used: 0.95 

$contrasts
Site = KL:
 contrast    estimate     SE df t.ratio p.value
 AC10 - AC12   0.1069 0.0381 18   2.804  0.0301
 AC10 - AC8    0.1952 0.0381 18   5.122  0.0002
 AC12 - AC8    0.0883 0.0381 18   2.318  0.0788

Site = SS:
 contrast    estimate     SE df t.ratio p.value
 AC10 - AC12   0.0157 0.0381 18   0.411  0.9115
 AC10 - AC8    0.1572 0.0381 18   4.124  0.0018
 AC12 - AC8    0.1415 0.0381 18   3.713  0.0043

P value adjustment: tukey method for comparing a family of 3 estimates 

#Sites within Genotypes
emmeans(Tol_ScoreF_M4_lm, pairwise~Site | Genotype)

$emmeans
Genotype = AC10:
 Site emmean    SE df lower.CL upper.CL
 KL    0.982 0.027 18    0.926    1.039
 SS    0.994 0.027 18    0.937    1.051

Genotype = AC12:
 Site emmean    SE df lower.CL upper.CL
 KL    0.876 0.027 18    0.819    0.932
 SS    0.978 0.027 18    0.922    1.035

Genotype = AC8:
 Site emmean    SE df lower.CL upper.CL
 KL    0.787 0.027 18    0.731    0.844
 SS    0.837 0.027 18    0.780    0.893

Confidence level used: 0.95 

$contrasts
Genotype = AC10:
 contrast estimate     SE df t.ratio p.value
 KL - SS   -0.0115 0.0381 18  -0.302  0.7663

Genotype = AC12:
 contrast estimate     SE df t.ratio p.value
 KL - SS   -0.1027 0.0381 18  -2.694  0.0148

Genotype = AC8:
 contrast estimate     SE df t.ratio p.value
 KL - SS   -0.0495 0.0381 18  -1.299  0.2103

##Save p-values

#Genotypes within Sites
Tol_ScoreF_M4_lm.geno<-data.frame(emmeans(Tol_ScoreF_M4_lm, pairwise~Genotype | Site)$contrasts)
Tol_ScoreF_M4_lm.geno<-Tol_ScoreF_M4_lm.geno %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_ScoreF_M4_lm.geno$group1<-paste(Tol_ScoreF_M4_lm.geno$Site, Tol_ScoreF_M4_lm.geno$group1, sep="_")
Tol_ScoreF_M4_lm.geno$group2<-paste(Tol_ScoreF_M4_lm.geno$Site, Tol_ScoreF_M4_lm.geno$group2, sep="_")

#Sites within Genotypes
Tol_ScoreF_M4_lm.site<-data.frame(emmeans(Tol_ScoreF_M4_lm, pairwise~Site | Genotype)$contrasts)
Tol_ScoreF_M4_lm.site<-Tol_ScoreF_M4_lm.site %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_ScoreF_M4_lm.site$group1<-paste(Tol_ScoreF_M4_lm.site$group1, Tol_ScoreF_M4_lm.site$Genotype, sep="_")
Tol_ScoreF_M4_lm.site$group2<-paste(Tol_ScoreF_M4_lm.site$group2, Tol_ScoreF_M4_lm.site$Genotype, sep="_")

#Full list of p-values
Tol_ScoreF_M4_lm.p<-rbind(Tol_ScoreF_M4_lm.geno[,c(1:2,4:8)], Tol_ScoreF_M4_lm.site[,c(1:2,4:8)])
Tol_ScoreF_M4_lm.p<-Tol_ScoreF_M4_lm.p %>% dplyr::rename( p = p.value)

#Add Significance Levels
Tol_ScoreF_M4_lm.p$Sig<-ifelse(Tol_ScoreF_M4_lm.p$p<0.001, "***", ifelse(Tol_ScoreF_M4_lm.p$p<0.01, "**", ifelse(Tol_ScoreF_M4_lm.p$p<0.05, "*", NA)))

#Specify Response and Timepoint
Tol_ScoreF_M4_lm.p$Response<-rep("Color_FullSet", nrow(Tol_ScoreF_M4_lm.p))
Tol_ScoreF_M4_lm.p$TimeP<-rep("M4", nrow(Tol_ScoreF_M4_lm.p))

Plot Retention by Site and Genotype

##Summary statistics by Site and Genotype
Tol_ScoreF_M4_SG<-summarySE(TolData_M4, measurevar="Score_Full.prop", groupvars=c("Site.Geno", "Site", "Genotype"), na.rm=TRUE)

##Plot Average Retention across Treatments
Tol_ScoreF_M4_SG.plot<-ggplot(Tol_ScoreF_M4_SG, aes(x=Site.Geno, y=Score_Full.prop, colour=Genotype)) + 
  geom_errorbar(aes(ymin=Score_Full.prop-se, ymax=Score_Full.prop+se), width=cap.sz, linewidth=bar.sz)+
  geom_point(size=point.sz)+ 
   ggtitle("Color Retention Full Set")+
  theme_classic()+
  theme( axis.title.x = element_text(size = axis.title.sz), 
         axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"),
        axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz),
        legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"),
        legend.position="bottom", 
        legend.direction="horizontal",
        plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
  labs(x="", y="Proportion Retained")+
  ylim(0, 1.5)+ 
  scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
  scale_color_manual(values = Geno.colors.o)+ 
  stat_pvalue_manual(data=Tol_ScoreF_M4_lm.p,  y.position=1.1, step.increase=0.22, label="Sig", hide.ns=TRUE); Tol_ScoreF_M4_SG.plot

Thermal Tolerance Color by Timepoint

Full Set

Run Model

##Check normality
hist(TolData$Score_TP.prop)

shapiro.test(TolData$Score_TP.prop)

    Shapiro-Wilk normality test

data:  TolData$Score_TP.prop
W = 0.9161, p-value = 0.0001986

#Not Normal

##Try square transformation
hist((TolData$Score_TP.prop)^2)

shapiro.test((TolData$Score_TP.prop)^2)

    Shapiro-Wilk normality test

data:  (TolData$Score_TP.prop)^2
W = 0.92506, p-value = 0.0004895

#Not normal 

##Try cubed transformation
hist((TolData$Score_TP.prop)^3)

shapiro.test((TolData$Score_TP.prop)^3)

    Shapiro-Wilk normality test

data:  (TolData$Score_TP.prop)^3
W = 0.92936, p-value = 0.0007674

#Not normal 


##Model as a function of Site and Genotype and Timepoint
##Model with no transformation and check residuals
Tol_ScoreTP_lm<-lm(Score_TP.prop~Site+Genotype+TimeP+ Site:Genotype + Site:TimeP + Genotype:TimeP, data=TolData)

Check Residuals

##Check Normality of Residuals
#Distribution 
plot(density(resid(Tol_ScoreTP_lm)))

#Q-Q plot
qqnorm(resid(Tol_ScoreTP_lm)); qqline(resid(Tol_ScoreTP_lm))

##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_ScoreTP_lm), resid(Tol_ScoreTP_lm))

Residuals are not great, need to check for other modeling options for Color by Timepoint.

Model Results

#Model Results
summary(Tol_ScoreTP_lm)

Call:
lm(formula = Score_TP.prop ~ Site + Genotype + TimeP + Site:Genotype + 
    Site:TimeP + Genotype:TimeP, data = TolData)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.133294 -0.019060  0.000725  0.026225  0.122906 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)         0.883494   0.006319 139.809  < 2e-16 ***
Site.L              0.038546   0.008912   4.325 6.47e-05 ***
Genotype.L         -0.102221   0.010797  -9.468 3.81e-13 ***
Genotype.Q          0.004971   0.011092   0.448 0.655784    
TimeP.L             0.068495   0.010797   6.344 4.45e-08 ***
TimeP.Q             0.021265   0.011092   1.917 0.060411 .  
Site.L:Genotype.L   0.014004   0.015269   0.917 0.363050    
Site.L:Genotype.Q  -0.008574   0.015625  -0.549 0.585424    
Site.L:TimeP.L     -0.036792   0.015269  -2.410 0.019344 *  
Site.L:TimeP.Q     -0.004886   0.015625  -0.313 0.755697    
Genotype.L:TimeP.L  0.029281   0.018818   1.556 0.125448    
Genotype.Q:TimeP.L  0.024642   0.018582   1.326 0.190292    
Genotype.L:TimeP.Q  0.015340   0.018582   0.826 0.412653    
Genotype.Q:TimeP.Q -0.072757   0.019821  -3.671 0.000548 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.05222 on 55 degrees of freedom
Multiple R-squared:  0.763, Adjusted R-squared:  0.7069 
F-statistic: 13.62 on 13 and 55 DF,  p-value: 8.098e-13

anova(Tol_ScoreTP_lm)

Analysis of Variance Table

Response: Score_TP.prop
               Df   Sum Sq  Mean Sq F value    Pr(>F)    
Site            1 0.057199 0.057199 20.9741 2.709e-05 ***
Genotype        2 0.238412 0.119206 43.7110 4.331e-12 ***
TimeP           2 0.117103 0.058552 21.4700 1.284e-07 ***
Site:Genotype   2 0.003221 0.001611  0.5906  0.557481    
Site:TimeP      2 0.017283 0.008642  3.1687  0.049832 *  
Genotype:TimeP  4 0.049580 0.012395  4.5450  0.003041 ** 
Residuals      55 0.149993 0.002727                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#Effect Size of Predictors
eta_squared(Tol_ScoreTP_lm, partial=FALSE)

# Effect Size for ANOVA (Type I)

Parameter      |     Eta2 |       95% CI
----------------------------------------
Site           |     0.09 | [0.01, 1.00]
Genotype       |     0.38 | [0.20, 1.00]
TimeP          |     0.19 | [0.04, 1.00]
Site:Genotype  | 5.09e-03 | [0.00, 1.00]
Site:TimeP     |     0.03 | [0.00, 1.00]
Genotype:TimeP |     0.08 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
Tol_ScoreTP_lm.res<-data.frame(anova(Tol_ScoreTP_lm))
Tol_ScoreTP_lm.res$Predictor<-rownames(Tol_ScoreTP_lm.res)
Tol_ScoreTP_lm.res$EtaSq<-c(eta_squared(Tol_ScoreTP_lm, partial=FALSE)$Eta2, NA)
Tol_ScoreTP_lm.res$Response<-rep("Color_TP", nrow(Tol_ScoreTP_lm.res))
Tol_ScoreTP_lm.res<-Tol_ScoreTP_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")

Strong influence of Timepoint, including interactions with variables of interest, Site and Genotype. Will analyze each Timepoint individually.

Color Score by Timepoint Summer 1 Timepoint W2

Run Model

##Check normality
hist(TolData_W2$Score_TP.prop)

shapiro.test(TolData_W2$Score_TP.prop)

    Shapiro-Wilk normality test

data:  TolData_W2$Score_TP.prop
W = 0.95832, p-value = 0.4302

#Normal

##Model as a function of Site and Genotype
Tol_ScoreTP_W2_lm<-lm(Score_TP.prop~Site+Genotype+Site:Genotype, data=TolData_W2)

Check Residuals

##Check Normality of Residuals
#Distribution 
plot(density(resid(Tol_ScoreTP_W2_lm)))

#Q-Q plot
qqnorm(resid(Tol_ScoreTP_W2_lm)); qqline(resid(Tol_ScoreTP_W2_lm))

##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_ScoreTP_W2_lm), resid(Tol_ScoreTP_W2_lm))

Model Results

Overall

#Model Results
summary(Tol_ScoreTP_W2_lm)

Call:
lm(formula = Score_TP.prop ~ Site + Genotype + Site:Genotype, 
    data = TolData_W2)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.117075 -0.028442 -0.002075  0.023658  0.099425 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)        0.843522   0.013080  64.488  < 2e-16 ***
Site.L             0.062257   0.018498   3.366  0.00367 ** 
Genotype.L        -0.117129   0.022952  -5.103 8.83e-05 ***
Genotype.Q        -0.042426   0.022356  -1.898  0.07485 .  
Site.L:Genotype.L -0.002071   0.032459  -0.064  0.94987    
Site.L:Genotype.Q  0.004044   0.031616   0.128  0.89972    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.06237 on 17 degrees of freedom
Multiple R-squared:  0.7169,    Adjusted R-squared:  0.6336 
F-statistic:  8.61 on 5 and 17 DF,  p-value: 0.0003242

anova(Tol_ScoreTP_W2_lm)

Analysis of Variance Table

Response: Score_TP.prop
              Df   Sum Sq  Mean Sq F value    Pr(>F)    
Site           1 0.054151 0.054151 13.9200 0.0016623 ** 
Genotype       2 0.113242 0.056621 14.5549 0.0002073 ***
Site:Genotype  2 0.000083 0.000041  0.0106 0.9894534    
Residuals     17 0.066133 0.003890                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#Effect Size of Predictors
eta_squared(Tol_ScoreTP_W2_lm, partial=FALSE)

# Effect Size for ANOVA (Type I)

Parameter     |     Eta2 |       95% CI
---------------------------------------
Site          |     0.23 | [0.01, 1.00]
Genotype      |     0.48 | [0.15, 1.00]
Site:Genotype | 3.53e-04 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
Tol_ScoreTP_W2_lm.res<-data.frame(anova(Tol_ScoreTP_W2_lm))
Tol_ScoreTP_W2_lm.res$Predictor<-rownames(Tol_ScoreTP_W2_lm.res)
Tol_ScoreTP_W2_lm.res$EtaSq<-c(eta_squared(Tol_ScoreTP_W2_lm, partial=FALSE)$Eta2, NA)
Tol_ScoreTP_W2_lm.res$Response<-rep("Color_TP", nrow(Tol_ScoreTP_W2_lm.res))
Tol_ScoreTP_W2_lm.res$TimeP<-rep("W2", nrow(Tol_ScoreTP_W2_lm.res))
Tol_ScoreTP_W2_lm.res<-Tol_ScoreTP_W2_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")

Significant effect of Site and Genotype. Still checking Site*Genotype for comparability across Timepoints.

Pairwise

#Pairwise comparisons across:

#Genotypes within Sites
emmeans(Tol_ScoreTP_W2_lm, pairwise~Genotype | Site)

$emmeans
Site = KL:
 Genotype emmean     SE df lower.CL upper.CL
 AC10      0.863 0.0312 17    0.797    0.929
 AC12      0.836 0.0312 17    0.771    0.902
 AC8       0.699 0.0312 17    0.633    0.765

Site = SS:
 Genotype emmean     SE df lower.CL upper.CL
 AC10      0.955 0.0312 17    0.889    1.021
 AC12      0.920 0.0312 17    0.854    0.986
 AC8       0.788 0.0360 17    0.712    0.864

Confidence level used: 0.95 

$contrasts
Site = KL:
 contrast    estimate     SE df t.ratio p.value
 AC10 - AC12   0.0263 0.0441 17   0.597  0.8236
 AC10 - AC8    0.1636 0.0441 17   3.709  0.0047
 AC12 - AC8    0.1373 0.0441 17   3.112  0.0166

Site = SS:
 contrast    estimate     SE df t.ratio p.value
 AC10 - AC12   0.0354 0.0441 17   0.803  0.7065
 AC10 - AC8    0.1677 0.0476 17   3.521  0.0070
 AC12 - AC8    0.1323 0.0476 17   2.778  0.0328

P value adjustment: tukey method for comparing a family of 3 estimates 

#Sites within Genotypes
emmeans(Tol_ScoreTP_W2_lm, pairwise~Site | Genotype)

$emmeans
Genotype = AC10:
 Site emmean     SE df lower.CL upper.CL
 KL    0.863 0.0312 17    0.797    0.929
 SS    0.955 0.0312 17    0.889    1.021

Genotype = AC12:
 Site emmean     SE df lower.CL upper.CL
 KL    0.836 0.0312 17    0.771    0.902
 SS    0.920 0.0312 17    0.854    0.986

Genotype = AC8:
 Site emmean     SE df lower.CL upper.CL
 KL    0.699 0.0312 17    0.633    0.765
 SS    0.788 0.0360 17    0.712    0.864

Confidence level used: 0.95 

$contrasts
Genotype = AC10:
 contrast estimate     SE df t.ratio p.value
 KL - SS   -0.0925 0.0441 17  -2.096  0.0513

Genotype = AC12:
 contrast estimate     SE df t.ratio p.value
 KL - SS   -0.0834 0.0441 17  -1.890  0.0759

Genotype = AC8:
 contrast estimate     SE df t.ratio p.value
 KL - SS   -0.0883 0.0476 17  -1.854  0.0812

##Save p-values

#Genotypes within Sites
Tol_ScoreTP_W2_lm.geno<-data.frame(emmeans(Tol_ScoreTP_W2_lm, pairwise~Genotype | Site)$contrasts)
Tol_ScoreTP_W2_lm.geno<-Tol_ScoreTP_W2_lm.geno %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_ScoreTP_W2_lm.geno$group1<-paste(Tol_ScoreTP_W2_lm.geno$Site, Tol_ScoreTP_W2_lm.geno$group1, sep="_")
Tol_ScoreTP_W2_lm.geno$group2<-paste(Tol_ScoreTP_W2_lm.geno$Site, Tol_ScoreTP_W2_lm.geno$group2, sep="_")

#Sites within Genotypes
Tol_ScoreTP_W2_lm.site<-data.frame(emmeans(Tol_ScoreTP_W2_lm, pairwise~Site | Genotype)$contrasts)
Tol_ScoreTP_W2_lm.site<-Tol_ScoreTP_W2_lm.site %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_ScoreTP_W2_lm.site$group1<-paste(Tol_ScoreTP_W2_lm.site$group1, Tol_ScoreTP_W2_lm.site$Genotype, sep="_")
Tol_ScoreTP_W2_lm.site$group2<-paste(Tol_ScoreTP_W2_lm.site$group2, Tol_ScoreTP_W2_lm.site$Genotype, sep="_")

#Full list of p-values
Tol_ScoreTP_W2_lm.p<-rbind(Tol_ScoreTP_W2_lm.geno[,c(1:2,4:8)], Tol_ScoreTP_W2_lm.site[,c(1:2,4:8)])
Tol_ScoreTP_W2_lm.p<-Tol_ScoreTP_W2_lm.p %>% dplyr::rename( p = p.value)

#Add Significance Levels
Tol_ScoreTP_W2_lm.p$Sig<-ifelse(Tol_ScoreTP_W2_lm.p$p<0.001, "***", ifelse(Tol_ScoreTP_W2_lm.p$p<0.01, "**", ifelse(Tol_ScoreTP_W2_lm.p$p<0.05, "*", NA)))

#Specify Response and Timepoint
Tol_ScoreTP_W2_lm.p$Response<-rep("Color_TP", nrow(Tol_ScoreTP_W2_lm.p))
Tol_ScoreTP_W2_lm.p$TimeP<-rep("W2", nrow(Tol_ScoreTP_W2_lm.p))

Plot Retention by Site and Genotype

##Summary statistics by Site and Genotype
Tol_ScoreTP_W2_SG<-summarySE(TolData_W2, measurevar="Score_TP.prop", groupvars=c("Site.Geno", "Site", "Genotype"), na.rm=TRUE)

##Plot Average Retention across Treatments
Tol_ScoreTP_W2_SG.plot<-ggplot(Tol_ScoreTP_W2_SG, aes(x=Site.Geno, y=Score_TP.prop, colour=Genotype)) + 
  geom_errorbar(aes(ymin=Score_TP.prop-se, ymax=Score_TP.prop+se), width=cap.sz, linewidth=bar.sz)+
  geom_point(size=point.sz)+ 
   ggtitle("Color Retention by Timepoint")+
  theme_classic()+
  theme( axis.title.x = element_text(size = axis.title.sz), 
         axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"),
        axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz),
        legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"),
        legend.position="bottom", 
        legend.direction="horizontal",
        plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
  labs(x="", y="Proportion Retained")+
  ylim(0, 1.5)+ 
  scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
  scale_color_manual(values = Geno.colors.o)+ 
  stat_pvalue_manual(data=Tol_ScoreTP_W2_lm.p,  y.position=1, step.increase=0.45, label="Sig", hide.ns=TRUE); Tol_ScoreTP_W2_SG.plot

Color Score by Timepoint Summer 2 Timepoint M1

Run Model

##Check normality
hist(TolData_M1$Score_TP.prop)

shapiro.test(TolData_M1$Score_TP.prop)

    Shapiro-Wilk normality test

data:  TolData_M1$Score_TP.prop
W = 0.9028, p-value = 0.03383

#Not normal

##Try square transformation
hist((TolData_M1$Score_TP.prop)^2)

shapiro.test((TolData_M1$Score_TP.prop)^2)

    Shapiro-Wilk normality test

data:  (TolData_M1$Score_TP.prop)^2
W = 0.90185, p-value = 0.03238

#Not normal

##Try cubed transformation
hist((TolData_M1$Score_TP.prop)^3)

shapiro.test((TolData_M1$Score_TP.prop)^3)

    Shapiro-Wilk normality test

data:  (TolData_M1$Score_TP.prop)^3
W = 0.89818, p-value = 0.02735

#Not normal

##Model as a function of Site and Genotype
##Model with no transformation and check residuals
Tol_ScoreTP_M1_lm<-lm(Score_TP.prop~Site+Genotype+Site:Genotype, data=TolData_M1)

Check Residuals

##Check Normality of Residuals
#Distribution 
plot(density(resid(Tol_ScoreTP_M1_lm)))

#Q-Q plot
qqnorm(resid(Tol_ScoreTP_M1_lm)); qqline(resid(Tol_ScoreTP_M1_lm))

##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_ScoreTP_M1_lm), resid(Tol_ScoreTP_M1_lm))

Model Results

Overall

#Model Results
summary(Tol_ScoreTP_M1_lm)

Call:
lm(formula = Score_TP.prop ~ Site + Genotype + Site:Genotype, 
    data = TolData_M1)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.125200 -0.020725 -0.001675  0.015662  0.131000 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)        0.866131   0.012681  68.303  < 2e-16 ***
Site.L             0.041582   0.017933   2.319    0.034 *  
Genotype.L        -0.114746   0.020837  -5.507 4.78e-05 ***
Genotype.Q         0.064377   0.023036   2.795    0.013 *  
Site.L:Genotype.L  0.050375   0.029467   1.710    0.107    
Site.L:Genotype.Q  0.004277   0.032577   0.131    0.897    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.05893 on 16 degrees of freedom
Multiple R-squared:  0.7446,    Adjusted R-squared:  0.6647 
F-statistic: 9.327 on 5 and 16 DF,  p-value: 0.0002606

anova(Tol_ScoreTP_M1_lm)

Analysis of Variance Table

Response: Score_TP.prop
              Df   Sum Sq  Mean Sq F value    Pr(>F)    
Site           1 0.019311 0.019311  5.5598   0.03144 *  
Genotype       2 0.132460 0.066230 19.0683 5.822e-05 ***
Site:Genotype  2 0.010210 0.005105  1.4698   0.25940    
Residuals     16 0.055573 0.003473                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#Effect Size of Predictors
eta_squared(Tol_ScoreTP_M1_lm, partial=FALSE)

# Effect Size for ANOVA (Type I)

Parameter     | Eta2 |       95% CI
-----------------------------------
Site          | 0.09 | [0.00, 1.00]
Genotype      | 0.61 | [0.30, 1.00]
Site:Genotype | 0.05 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
Tol_ScoreTP_M1_lm.res<-data.frame(anova(Tol_ScoreTP_M1_lm))
Tol_ScoreTP_M1_lm.res$Predictor<-rownames(Tol_ScoreTP_M1_lm.res)
Tol_ScoreTP_M1_lm.res$EtaSq<-c(eta_squared(Tol_ScoreTP_M1_lm, partial=FALSE)$Eta2, NA)
Tol_ScoreTP_M1_lm.res$Response<-rep("Color_TP", nrow(Tol_ScoreTP_M1_lm.res))
Tol_ScoreTP_M1_lm.res$TimeP<-rep("M1", nrow(Tol_ScoreTP_M1_lm.res))
Tol_ScoreTP_M1_lm.res<-Tol_ScoreTP_M1_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")

Significant effect of Site and Genotype. Still checking Site*Genotype for comparability across Timepoints.

Pairwise

#Pairwise comparisons across:

#Genotypes within Sites
emmeans(Tol_ScoreTP_M1_lm, pairwise~Genotype | Site)

$emmeans
Site = KL:
 Genotype emmean     SE df lower.CL upper.CL
 AC10      0.968 0.0295 16    0.906    1.031
 AC12      0.787 0.0340 16    0.715    0.859
 AC8       0.755 0.0295 16    0.693    0.818

Site = SS:
 Genotype emmean     SE df lower.CL upper.CL
 AC10      0.979 0.0295 16    0.917    1.041
 AC12      0.841 0.0340 16    0.768    0.913
 AC8       0.867 0.0295 16    0.805    0.930

Confidence level used: 0.95 

$contrasts
Site = KL:
 contrast    estimate     SE df t.ratio p.value
 AC10 - AC12   0.1815 0.0450 16   4.032  0.0026
 AC10 - AC8    0.2127 0.0417 16   5.103  0.0003
 AC12 - AC8    0.0312 0.0450 16   0.693  0.7710

Site = SS:
 contrast    estimate     SE df t.ratio p.value
 AC10 - AC12   0.1385 0.0450 16   3.077  0.0187
 AC10 - AC8    0.1119 0.0417 16   2.685  0.0408
 AC12 - AC8   -0.0266 0.0450 16  -0.591  0.8269

P value adjustment: tukey method for comparing a family of 3 estimates 

#Sites within Genotypes
emmeans(Tol_ScoreTP_M1_lm, pairwise~Site | Genotype)

$emmeans
Genotype = AC10:
 Site emmean     SE df lower.CL upper.CL
 KL    0.968 0.0295 16    0.906    1.031
 SS    0.979 0.0295 16    0.917    1.041

Genotype = AC12:
 Site emmean     SE df lower.CL upper.CL
 KL    0.787 0.0340 16    0.715    0.859
 SS    0.841 0.0340 16    0.768    0.913

Genotype = AC8:
 Site emmean     SE df lower.CL upper.CL
 KL    0.755 0.0295 16    0.693    0.818
 SS    0.867 0.0295 16    0.805    0.930

Confidence level used: 0.95 

$contrasts
Genotype = AC10:
 contrast estimate     SE df t.ratio p.value
 KL - SS   -0.0109 0.0417 16  -0.262  0.7970

Genotype = AC12:
 contrast estimate     SE df t.ratio p.value
 KL - SS   -0.0539 0.0481 16  -1.119  0.2795

Genotype = AC8:
 contrast estimate     SE df t.ratio p.value
 KL - SS   -0.1116 0.0417 16  -2.679  0.0165

##Save p-values

#Genotypes within Sites
Tol_ScoreTP_M1_lm.geno<-data.frame(emmeans(Tol_ScoreTP_M1_lm, pairwise~Genotype | Site)$contrasts)
Tol_ScoreTP_M1_lm.geno<-Tol_ScoreTP_M1_lm.geno %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_ScoreTP_M1_lm.geno$group1<-paste(Tol_ScoreTP_M1_lm.geno$Site, Tol_ScoreTP_M1_lm.geno$group1, sep="_")
Tol_ScoreTP_M1_lm.geno$group2<-paste(Tol_ScoreTP_M1_lm.geno$Site, Tol_ScoreTP_M1_lm.geno$group2, sep="_")

#Sites within Genotypes
Tol_ScoreTP_M1_lm.site<-data.frame(emmeans(Tol_ScoreTP_M1_lm, pairwise~Site | Genotype)$contrasts)
Tol_ScoreTP_M1_lm.site<-Tol_ScoreTP_M1_lm.site %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_ScoreTP_M1_lm.site$group1<-paste(Tol_ScoreTP_M1_lm.site$group1, Tol_ScoreTP_M1_lm.site$Genotype, sep="_")
Tol_ScoreTP_M1_lm.site$group2<-paste(Tol_ScoreTP_M1_lm.site$group2, Tol_ScoreTP_M1_lm.site$Genotype, sep="_")

#Full list of p-values
Tol_ScoreTP_M1_lm.p<-rbind(Tol_ScoreTP_M1_lm.geno[,c(1:2,4:8)], Tol_ScoreTP_M1_lm.site[,c(1:2,4:8)])
Tol_ScoreTP_M1_lm.p<-Tol_ScoreTP_M1_lm.p %>% dplyr::rename( p = p.value)

#Add Significance Levels
Tol_ScoreTP_M1_lm.p$Sig<-ifelse(Tol_ScoreTP_M1_lm.p$p<0.001, "***", ifelse(Tol_ScoreTP_M1_lm.p$p<0.01, "**", ifelse(Tol_ScoreTP_M1_lm.p$p<0.05, "*", NA)))

#Specify Response and Timepoint
Tol_ScoreTP_M1_lm.p$Response<-rep("Color_TP", nrow(Tol_ScoreTP_M1_lm.p))
Tol_ScoreTP_M1_lm.p$TimeP<-rep("M1", nrow(Tol_ScoreTP_M1_lm.p))

Plot Retention by Site and Genotype

##Summary statistics by Site and Genotype
Tol_ScoreTP_M1_SG<-summarySE(TolData_M1, measurevar="Score_TP.prop", groupvars=c("Site.Geno", "Site", "Genotype"), na.rm=TRUE)

##Plot Average Retention across Treatments
Tol_ScoreTP_M1_SG.plot<-ggplot(Tol_ScoreTP_M1_SG, aes(x=Site.Geno, y=Score_TP.prop, colour=Genotype)) + 
  geom_errorbar(aes(ymin=Score_TP.prop-se, ymax=Score_TP.prop+se), width=cap.sz, linewidth=bar.sz)+
  geom_point(size=point.sz)+ 
   ggtitle("Color Retention by Timepoint")+
  theme_classic()+
  theme( axis.title.x = element_text(size = axis.title.sz), 
         axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"),
        axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz),
        legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"),
        legend.position="bottom", 
        legend.direction="horizontal",
        plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
  labs(x="", y="Proportion Retained")+
 ylim(0, 1.5)+ 
  scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
  scale_color_manual(values = Geno.colors.o)+ 
  stat_pvalue_manual(data=Tol_ScoreTP_M1_lm.p,  y.position=1.05, step.increase=0.3, label="Sig", hide.ns=TRUE); Tol_ScoreTP_M1_SG.plot

Color Score by Timepoint Winter Timepoint M4

Run Model

##Check normality
hist(TolData_M4$Score_TP.prop)

shapiro.test(TolData_M4$Score_TP.prop)

    Shapiro-Wilk normality test

data:  TolData_M4$Score_TP.prop
W = 0.85397, p-value = 0.002593

#Not Normal

##Try square transformation
hist((TolData_M4$Score_TP.prop)^2)

shapiro.test((TolData_M4$Score_TP.prop)^2)

    Shapiro-Wilk normality test

data:  (TolData_M4$Score_TP.prop)^2
W = 0.85439, p-value = 0.00264

#Not Normal

##Try cubed transformation
hist((TolData_M4$Score_TP.prop)^3)

shapiro.test((TolData_M4$Score_TP.prop)^3)

    Shapiro-Wilk normality test

data:  (TolData_M4$Score_TP.prop)^3
W = 0.85457, p-value = 0.002661

#Not Normal

##Model as a function of Site and Genotype
##Model with no transformation and check residuals
Tol_ScoreTP_M4_lm<-lm(Score_TP.prop~Site+Genotype+Site:Genotype, data=TolData_M4)

Check Residuals

##Check Normality of Residuals
#Distribution 
plot(density(resid(Tol_ScoreTP_M4_lm)))

#Q-Q plot
qqnorm(resid(Tol_ScoreTP_M4_lm)); qqline(resid(Tol_ScoreTP_M4_lm))

##Check Variance of Residuals across Fitted Values
plot(fitted(Tol_ScoreTP_M4_lm), resid(Tol_ScoreTP_M4_lm))

Model Results

Overall

#Model Results
summary(Tol_ScoreTP_M4_lm)

Call:
lm(formula = Score_TP.prop ~ Site + Genotype + Site:Genotype, 
    data = TolData_M4)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.070400 -0.009438  0.004100  0.011281  0.057100 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)        0.940608   0.006243 150.654  < 2e-16 ***
Site.L             0.010536   0.008830   1.193   0.2483    
Genotype.L        -0.075254   0.010814  -6.959 1.68e-06 ***
Genotype.Q        -0.007308   0.010814  -0.676   0.5078    
Site.L:Genotype.L -0.006950   0.015293  -0.454   0.6549    
Site.L:Genotype.Q -0.032086   0.015293  -2.098   0.0503 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.03059 on 18 degrees of freedom
Multiple R-squared:  0.7531,    Adjusted R-squared:  0.6846 
F-statistic: 10.98 on 5 and 18 DF,  p-value: 5.763e-05

anova(Tol_ScoreTP_M4_lm)

Analysis of Variance Table

Response: Score_TP.prop
              Df   Sum Sq   Mean Sq F value    Pr(>F)    
Site           1 0.001332 0.0013321  1.4238    0.2483    
Genotype       2 0.045732 0.0228662 24.4414 7.407e-06 ***
Site:Genotype  2 0.004311 0.0021557  2.3042    0.1285    
Residuals     18 0.016840 0.0009355                      
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#Effect Size of Predictors
eta_squared(Tol_ScoreTP_M4_lm, partial=FALSE)

# Effect Size for ANOVA (Type I)

Parameter     | Eta2 |       95% CI
-----------------------------------
Site          | 0.02 | [0.00, 1.00]
Genotype      | 0.67 | [0.41, 1.00]
Site:Genotype | 0.06 | [0.00, 1.00]

- One-sided CIs: upper bound fixed at [1.00].

##Save model results
Tol_ScoreTP_M4_lm.res<-data.frame(anova(Tol_ScoreTP_M4_lm))
Tol_ScoreTP_M4_lm.res$Predictor<-rownames(Tol_ScoreTP_M4_lm.res)
Tol_ScoreTP_M4_lm.res$EtaSq<-c(eta_squared(Tol_ScoreTP_M4_lm, partial=FALSE)$Eta2, NA)
Tol_ScoreTP_M4_lm.res$Response<-rep("Color_TP", nrow(Tol_ScoreTP_M4_lm.res))
Tol_ScoreTP_M4_lm.res$TimeP<-rep("M4", nrow(Tol_ScoreTP_M4_lm.res))
Tol_ScoreTP_M4_lm.res<-Tol_ScoreTP_M4_lm.res %>% dplyr::rename( p.value = "Pr..F.", DF= "Df")

Significant effect of Genotype. Still checking Site*Genotype for comparability across Timepoints.

Pairwise

#Pairwise comparisons across:

#Genotypes within Sites
emmeans(Tol_ScoreTP_M4_lm, pairwise~Genotype | Site)

$emmeans
Site = KL:
 Genotype emmean     SE df lower.CL upper.CL
 AC10      0.989 0.0153 18    0.957    1.021
 AC12      0.921 0.0153 18    0.888    0.953
 AC8       0.890 0.0153 18    0.858    0.922

Site = SS:
 Genotype emmean     SE df lower.CL upper.CL
 AC10      0.993 0.0153 18    0.960    1.025
 AC12      0.973 0.0153 18    0.940    1.005
 AC8       0.879 0.0153 18    0.847    0.911

Confidence level used: 0.95 

$contrasts
Site = KL:
 contrast    estimate     SE df t.ratio p.value
 AC10 - AC12   0.0686 0.0216 18   3.171  0.0139
 AC10 - AC8    0.0995 0.0216 18   4.599  0.0006
 AC12 - AC8    0.0309 0.0216 18   1.429  0.3478

Site = SS:
 contrast    estimate     SE df t.ratio p.value
 AC10 - AC12   0.0199 0.0216 18   0.922  0.6335
 AC10 - AC8    0.1134 0.0216 18   5.242  0.0002
 AC12 - AC8    0.0934 0.0216 18   4.320  0.0011

P value adjustment: tukey method for comparing a family of 3 estimates 

#Sites within Genotypes
emmeans(Tol_ScoreTP_M4_lm, pairwise~Site | Genotype)

$emmeans
Genotype = AC10:
 Site emmean     SE df lower.CL upper.CL
 KL    0.989 0.0153 18    0.957    1.021
 SS    0.993 0.0153 18    0.960    1.025

Genotype = AC12:
 Site emmean     SE df lower.CL upper.CL
 KL    0.921 0.0153 18    0.888    0.953
 SS    0.973 0.0153 18    0.940    1.005

Genotype = AC8:
 Site emmean     SE df lower.CL upper.CL
 KL    0.890 0.0153 18    0.858    0.922
 SS    0.879 0.0153 18    0.847    0.911

Confidence level used: 0.95 

$contrasts
Genotype = AC10:
 contrast estimate     SE df t.ratio p.value
 KL - SS  -0.00332 0.0216 18  -0.154  0.8795

Genotype = AC12:
 contrast estimate     SE df t.ratio p.value
 KL - SS  -0.05195 0.0216 18  -2.402  0.0273

Genotype = AC8:
 contrast estimate     SE df t.ratio p.value
 KL - SS   0.01057 0.0216 18   0.489  0.6308

##Save p-values

#Genotypes within Sites
Tol_ScoreTP_M4_lm.geno<-data.frame(emmeans(Tol_ScoreTP_M4_lm, pairwise~Genotype | Site)$contrasts)
Tol_ScoreTP_M4_lm.geno<-Tol_ScoreTP_M4_lm.geno %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_ScoreTP_M4_lm.geno$group1<-paste(Tol_ScoreTP_M4_lm.geno$Site, Tol_ScoreTP_M4_lm.geno$group1, sep="_")
Tol_ScoreTP_M4_lm.geno$group2<-paste(Tol_ScoreTP_M4_lm.geno$Site, Tol_ScoreTP_M4_lm.geno$group2, sep="_")

#Sites within Genotypes
Tol_ScoreTP_M4_lm.site<-data.frame(emmeans(Tol_ScoreTP_M4_lm, pairwise~Site | Genotype)$contrasts)
Tol_ScoreTP_M4_lm.site<-Tol_ScoreTP_M4_lm.site %>% separate(col=contrast, into=c("group1", "group2"), sep=" - ", remove=TRUE)
Tol_ScoreTP_M4_lm.site$group1<-paste(Tol_ScoreTP_M4_lm.site$group1, Tol_ScoreTP_M4_lm.site$Genotype, sep="_")
Tol_ScoreTP_M4_lm.site$group2<-paste(Tol_ScoreTP_M4_lm.site$group2, Tol_ScoreTP_M4_lm.site$Genotype, sep="_")

#Full list of p-values
Tol_ScoreTP_M4_lm.p<-rbind(Tol_ScoreTP_M4_lm.geno[,c(1:2,4:8)], Tol_ScoreTP_M4_lm.site[,c(1:2,4:8)])
Tol_ScoreTP_M4_lm.p<-Tol_ScoreTP_M4_lm.p %>% dplyr::rename( p = p.value)

#Add Significance Levels
Tol_ScoreTP_M4_lm.p$Sig<-ifelse(Tol_ScoreTP_M4_lm.p$p<0.001, "***", ifelse(Tol_ScoreTP_M4_lm.p$p<0.01, "**", ifelse(Tol_ScoreTP_M4_lm.p$p<0.05, "*", NA)))

#Specify Response and Timepoint
Tol_ScoreTP_M4_lm.p$Response<-rep("Color_TP", nrow(Tol_ScoreTP_M4_lm.p))
Tol_ScoreTP_M4_lm.p$TimeP<-rep("M4", nrow(Tol_ScoreTP_M4_lm.p))

Plot Retention by Site and Genotype

##Summary statistics by Site and Genotype
Tol_ScoreTP_M4_SG<-summarySE(TolData_M4, measurevar="Score_TP.prop", groupvars=c("Site.Geno", "Site", "Genotype"), na.rm=TRUE)

##Plot Average Retention across Treatments
Tol_ScoreTP_M4_SG.plot<-ggplot(Tol_ScoreTP_M4_SG, aes(x=Site.Geno, y=Score_TP.prop, colour=Genotype)) + 
  geom_errorbar(aes(ymin=Score_TP.prop-se, ymax=Score_TP.prop+se), width=cap.sz, linewidth=bar.sz)+
  geom_point(size=point.sz)+ 
   ggtitle("Color Retention by Timepoint")+
  theme_classic()+
  theme( axis.title.x = element_text(size = axis.title.sz), 
         axis.title.y = element_text(size = axis.title.sz), 
        axis.text.x=element_text(size=axis.txt.sz, colour="black"),
        axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
        legend.text=element_text(size=leg.txt.sz),
        legend.title=element_text(size=leg.title.sz), 
        legend.box.background = element_rect(color = "black"),
        legend.position="bottom", 
        legend.direction="horizontal",
        plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
  labs(x="", y="Proportion Retained")+
  ylim(0, 1.5)+ 
  scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
  scale_color_manual(values = Geno.colors.o)+ 
  stat_pvalue_manual(data=Tol_ScoreTP_M4_lm.p,  y.position=1.1, step.increase=0.4, label="Sig", hide.ns=TRUE); Tol_ScoreTP_M4_SG.plot

Contrasts Across Metrics

Combine Contrast Results

Creating a heatmap to compare the direction and significance of pairwise comparisons across Tolerance metrics. Positive estimates indicate that Group 1 > Group 2. P values of less than 0.05 are considered significant.

##Combine Results
Tol_contrasts<-rbind(Tol_Chl_W2_lm.p, Tol_Chl_M1_lm.p, Tol_Chl_M4_lm.p, 
                     Tol_Sym_W2_lm.p, Tol_Sym_M1_lm.p, Tol_Sym_M4_lm.p, 
                     Tol_ScoreF_W2_lm.p, Tol_ScoreF_M1_lm.p, Tol_ScoreF_M4_lm.p, 
                     Tol_ScoreTP_W2_lm.p, Tol_ScoreTP_M1_lm.p, Tol_ScoreTP_M4_lm.p)

##Create Contrast Column
Tol_contrasts$Contrast<-paste(Tol_contrasts$group1, Tol_contrasts$group2, sep=" v. ")

##Add Set column with Contrast and Timepoint
Tol_contrasts$Set<-paste(Tol_contrasts$TimeP, Tol_contrasts$Contrast, sep=" ")

##Re-order by Seasonal Timepoint
Tol_contrasts$TimeP<-factor(Tol_contrasts$TimeP, levels=c("W2", "M1", "M4"), ordered=TRUE)
Tol_contrasts<- Tol_contrasts %>% arrange(desc(as.numeric(TimeP)))

Tol_contrasts$Set<-factor(Tol_contrasts$Set, levels=c(unique(Tol_contrasts$Set)), ordered=TRUE)

Standardize Response Scale

Convert Estimate to the Response scale (instead of log +1 scale) for models where the response was log transformed

Tol_contrasts$Estimate<-Tol_contrasts$estimate

#Chl M4
Tol_contrasts$Estimate[which(Tol_contrasts$Response=="Chlorophyll" & Tol_contrasts$TimeP=="M4")]<-c(exp(Tol_contrasts$estimate[which(Tol_contrasts$Response=="Chlorophyll" & Tol_contrasts$TimeP=="M4")])-1)

#Sym M1
Tol_contrasts$Estimate[which(Tol_contrasts$Response=="Symbionts" & Tol_contrasts$TimeP=="M1")]<-c(exp(Tol_contrasts$estimate[which(Tol_contrasts$Response=="Symbionts" & Tol_contrasts$TimeP=="M1")])-1)

Contrast Heatmaps

Summer 1 Timepoint W2

Tol_Pairs_W2.plot<-ggplot(data=Tol_contrasts[which(Tol_contrasts$TimeP=="W2"),], aes(x=Response, y=Contrast, fill=Estimate))+
  geom_tile()+
  geom_text(aes(label=Sig), size=sig.sz)+
  scale_fill_gradient2(low="#BA1E02FF", mid="white", high="#3BA0FDFF")+
  theme_bw()+
  ggtitle("Pairwise Contrasts Summer 1")+
  theme(axis.text.x=element_text(size=axis.title.sz, colour="black", angle=45, hjust=1),
        axis.text.y=element_text(size=axis.txt.sz-1, colour="black"),
        legend.text=element_text(size=leg.txt.sz),
        legend.title=element_text(size=leg.title.sz), 
        plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
  labs(x="", y="", fill="Estimate");Tol_Pairs_W2.plot

Summer 2 Timepoint M1

Tol_Pairs_M1.plot<-ggplot(data=Tol_contrasts[which(Tol_contrasts$TimeP=="M1"),], aes(x=Response, y=Contrast, fill=Estimate))+
  geom_tile()+
  geom_text(aes(label=Sig), size=sig.sz)+
  scale_fill_gradient2(low="#BA1E02FF", mid="white", high="#3BA0FDFF")+
  theme_bw()+
  ggtitle("Pairwise Contrasts Summer 2")+
  theme(axis.text.x=element_text(size=axis.title.sz, colour="black", angle=45, hjust=1),
        axis.text.y=element_text(size=axis.txt.sz-1, colour="black"),
        legend.text=element_text(size=leg.txt.sz),
        legend.title=element_text(size=leg.title.sz), 
        plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
  labs(x="", y="", fill="Estimate");Tol_Pairs_M1.plot

Winter Timepoint M4

Tol_Pairs_M4.plot<-ggplot(data=Tol_contrasts[which(Tol_contrasts$TimeP=="M4"),], aes(x=Response, y=Contrast, fill=Estimate))+
  geom_tile()+
  geom_text(aes(label=Sig), size=sig.sz)+
  scale_fill_gradient2(low="#BA1E02FF", mid="white", high="#3BA0FDFF")+
  theme_bw()+
  ggtitle("Pairwise Contrasts Month 4")+
  theme(axis.text.x=element_text(size=axis.title.sz, colour="black", angle=45, hjust=1),
        axis.text.y=element_text(size=axis.txt.sz-1, colour="black"),
        legend.text=element_text(size=leg.txt.sz),
        legend.title=element_text(size=leg.title.sz), 
        plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
  labs(x="", y="", fill="Estimate");Tol_Pairs_M4.plot

All Timepoints

Tol_Pairs.plot<-ggplot(data=Tol_contrasts, aes(x=Response, y=Set, fill=Estimate))+
  geom_tile()+
  geom_text(aes(label=Sig), size=sig.sz)+
  scale_fill_gradient2(low="#BA1E02FF", mid="white", high="#3BA0FDFF")+
  theme_bw()+
  ggtitle("Pairwise Contrasts")+
  theme(axis.text.x=element_text(size=axis.title.sz, colour="black", angle=45, hjust=1),
        axis.text.y=element_text(size=axis.txt.sz-1, colour="black"),
        legend.text=element_text(size=leg.txt.sz),
        legend.title=element_text(size=leg.title.sz), 
        plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
  labs(x="", y="", fill="Estimate");Tol_Pairs.plot

Significant Differences

Subsetting Contrasts to only include contrasts with significant differences in at least one of the thermal tolerance metrics (Chlorophyll or Symbiont or Color retention)

##Remove rows with non-significant p-values
Tol_contrasts_sig<-subset(Tol_contrasts, p < 0.05)

##Keep Sets where at least one metric shows a significant difference
Tol_Sets_sig<-data.frame(Set=c(unique(Tol_contrasts_sig$Set)))
Tol_contrasts_sig_Set<-merge(Tol_Sets_sig, Tol_contrasts, all.x=TRUE, all.y=FALSE)

Tol_Pairs_sig.plot<-ggplot(data=Tol_contrasts_sig_Set, aes(x=Response, y=Set, fill=Estimate))+
  geom_tile()+
  geom_text(aes(label=Sig), size=sig.sz)+
  scale_fill_gradient2(low="#BA1E02FF", mid="white", high="#3BA0FDFF")+
  theme_bw()+
  ggtitle("Significant Pairwise Contrasts")+
  theme(axis.text.x=element_text(size=axis.title.sz, colour="black", angle=45, hjust=1),
        axis.text.y=element_text(size=axis.txt.sz-1, colour="black"),
        legend.text=element_text(size=leg.txt.sz),
        legend.title=element_text(size=leg.title.sz), 
        plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
  labs(x="", y="", fill="Estimate");Tol_Pairs_sig.plot

Proportion of Agreement

Assessing the proportion of agreement between methods as the pairwise comparisons where two methods both find a significant difference out of the total pairwise comparisons. Calculating the proportion of agreement between Chlorophyll and Color Full Set, Chlorophyll and Color Timepoint, Symbionts and Color Full Set, Symbionts and Color Timepoint, and Chlorophyll and Symbionts. Comparing the proportion of agreement between methods of quantifying thermal tolerance.

names(Tol_contrasts)

 [1] "group1"   "group2"   "estimate" "SE"       "df"       "t.ratio"  "p"        "Sig"      "Response"
[10] "TimeP"    "Contrast" "Set"      "Estimate"

##Change long to short format
Tol_contrasts.agree<-Tol_contrasts %>% pivot_wider(names_from="Response", values_from=c("p", "Estimate"), id_cols=c("TimeP", "Contrast", "Set"))

##Convert p values to binary
#0: No significant difference (p>0.05) detected in that pairwise contrast
#1: Significant difference (p<0.05) detected in that pairwise contrast
Tol_contrasts.agree[,c(4:7)]<-ifelse(Tol_contrasts.agree[,c(4:7)]<0.05, 1, 0)

##Convert estimates to Greater or Less
#Greater: Positive estimates indicate A>B in the Contrast
#Less: Negative estimates indicate A<B in the Contrast
Tol_contrasts.agree[,c(8:11)]<-ifelse(Tol_contrasts.agree[,c(8:11)]<0, "Less", "Greater")

##Identify agreement
#1: Both metrics identified a significant difference in that pairwise contrast OR Neither metric identified a significant difference
#0: One metrics identified a significant difference in that pairwise contrast
#For each matching case (1), confirm matching Estimates
Tol_contrasts.agree$Chl.Sym<-ifelse(Tol_contrasts.agree$p_Chlorophyll==1 &
Tol_contrasts.agree$p_Chlorophyll==Tol_contrasts.agree$p_Symbionts &
Tol_contrasts.agree$Estimate_Chlorophyll==Tol_contrasts.agree$Estimate_Symbionts, 1, 
ifelse(Tol_contrasts.agree$p_Chlorophyll==0 & 
         Tol_contrasts.agree$p_Chlorophyll==Tol_contrasts.agree$p_Symbionts, 1, 0))

Tol_contrasts.agree$Chl.ColF<-ifelse(Tol_contrasts.agree$p_Chlorophyll==1 &                                     Tol_contrasts.agree$p_Chlorophyll==Tol_contrasts.agree$p_Color_FullSet &
Tol_contrasts.agree$Estimate_Chlorophyll==Tol_contrasts.agree$Estimate_Color_FullSet, 1, 
ifelse(Tol_contrasts.agree$p_Chlorophyll==0 &
         Tol_contrasts.agree$p_Chlorophyll==Tol_contrasts.agree$p_Color_FullSet, 1, 0))

Tol_contrasts.agree$Chl.ColTP<-ifelse(Tol_contrasts.agree$p_Chlorophyll==1 & 
Tol_contrasts.agree$p_Chlorophyll==Tol_contrasts.agree$p_Color_TP &                          Tol_contrasts.agree$Estimate_Chlorophyll==Tol_contrasts.agree$Estimate_Color_TP, 1, 
ifelse(Tol_contrasts.agree$p_Chlorophyll==0 &
         Tol_contrasts.agree$p_Chlorophyll==Tol_contrasts.agree$p_Color_TP, 1, 0))

Tol_contrasts.agree$Sym.ColF<-ifelse(Tol_contrasts.agree$p_Symbionts==1 &
Tol_contrasts.agree$p_Symbionts==Tol_contrasts.agree$p_Color_FullSet & Tol_contrasts.agree$Estimate_Symbionts==Tol_contrasts.agree$Estimate_Color_FullSet, 1,
ifelse(Tol_contrasts.agree$p_Symbionts==0 &
         Tol_contrasts.agree$p_Symbionts==Tol_contrasts.agree$p_Color_FullSet, 1, 0))

Tol_contrasts.agree$Sym.ColTP<-ifelse(Tol_contrasts.agree$p_Symbionts==1 &
Tol_contrasts.agree$p_Symbionts==Tol_contrasts.agree$p_Color_TP & Tol_contrasts.agree$Estimate_Symbionts==Tol_contrasts.agree$Estimate_Color_TP, 1,
ifelse(Tol_contrasts.agree$p_Symbionts==0 &
         Tol_contrasts.agree$p_Symbionts==Tol_contrasts.agree$p_Color_TP, 1, 0))

Tol_contrasts.agree$ColF.ColTP<-ifelse(Tol_contrasts.agree$p_Color_FullSet==1 &
Tol_contrasts.agree$p_Color_FullSet==Tol_contrasts.agree$p_Color_TP & Tol_contrasts.agree$Estimate_Color_FullSet==Tol_contrasts.agree$Estimate_Color_TP, 1, 
ifelse(Tol_contrasts.agree$p_Color_FullSet==0 &
         Tol_contrasts.agree$p_Color_FullSet==Tol_contrasts.agree$p_Color_TP, 1, 0))

##Convert to long to calculate instances of agreement
Tol_contrasts.agree.long<-Tol_contrasts.agree %>% pivot_longer(cols=c("Chl.Sym", "Chl.ColF", "Chl.ColTP", "Sym.ColF", "Sym.ColTP", "ColF.ColTP"), names_to="Compare", values_to="Agree")

##Sum Cases of Agreement
Tol_contrasts.agree.sum<-aggregate(Tol_contrasts.agree.long$Agree, list(Tol_contrasts.agree.long$Compare), sum)
names(Tol_contrasts.agree.sum)<-c("Compare", "Agree")

##Total Pairwise Comparisons
Tol_contrasts.agree.sum$Total<-aggregate(Tol_contrasts.agree.long$Agree, list(Tol_contrasts.agree.long$Compare), function(x) length(x))[,2]

Compare Proportion Agreement

prop.test(Tol_contrasts.agree.sum$Agree, Tol_contrasts.agree.sum$Total)

    6-sample test for equality of proportions without continuity correction

data:  Tol_contrasts.agree.sum$Agree out of Tol_contrasts.agree.sum$Total
X-squared = 6.0509, df = 5, p-value = 0.3013
alternative hypothesis: two.sided
sample estimates:
   prop 1    prop 2    prop 3    prop 4    prop 5    prop 6 
0.7777778 0.7777778 0.7037037 0.9259259 0.6666667 0.7407407 

No significant difference in the proportion of agreement in detecting pairwise differences in thermal tolerance across comparisons of traditional bleaching or color score methods (Pearson’s chi-squared test p-value = 0.3013)

Compare Agreement across Seasons

##Sum Cases of Agreement by Comparison and Timepoint
Tol_contrasts.agree.sum.tp<-aggregate(Tol_contrasts.agree.long$Agree, list(Tol_contrasts.agree.long$TimeP, Tol_contrasts.agree.long$Compare), sum)
names(Tol_contrasts.agree.sum.tp)<-c("TimeP", "Compare", "Agree")

##Total Pairwise Comparisons
Tol_contrasts.agree.sum.tp$Total<-aggregate(Tol_contrasts.agree.long$Agree, list(Tol_contrasts.agree.long$TimeP, Tol_contrasts.agree.long$Compare), function(x) length(x))[,3]

Chl vs Sym

##Subset by Comparison
Tol_contrasts.agree.Chl.Sym<-subset(Tol_contrasts.agree.sum.tp, Compare=="Chl.Sym")

##Compare Proportions between Timepoints
prop.test(Tol_contrasts.agree.Chl.Sym$Agree, Tol_contrasts.agree.Chl.Sym$Total)

Warning: Chi-squared approximation may be incorrect

    3-sample test for equality of proportions without continuity correction

data:  Tol_contrasts.agree.Chl.Sym$Agree out of Tol_contrasts.agree.Chl.Sym$Total
X-squared = 5.6842, df = 2, p-value = 0.0583
alternative hypothesis: two.sided
sample estimates:
   prop 1    prop 2    prop 3 
0.5555556 1.0000000 0.5555556 

Chl vs Color Full

##Subset by Comparison
Tol_contrasts.agree.Chl.ColF<-subset(Tol_contrasts.agree.sum.tp, Compare=="Chl.ColF")

##Compare Proportions between Timepoints
prop.test(Tol_contrasts.agree.Chl.ColF$Agree, Tol_contrasts.agree.Chl.ColF$Total)

Warning: Chi-squared approximation may be incorrect

    3-sample test for equality of proportions without continuity correction

data:  Tol_contrasts.agree.Chl.ColF$Agree out of Tol_contrasts.agree.Chl.ColF$Total
X-squared = 9, df = 2, p-value = 0.01111
alternative hypothesis: two.sided
sample estimates:
   prop 1    prop 2    prop 3 
1.0000000 0.8888889 0.4444444 

Chl vs Color Timepoint

##Subset by Comparison
Tol_contrasts.agree.Chl.ColTP<-subset(Tol_contrasts.agree.sum.tp, Compare=="Chl.ColTP")

##Compare Proportions between Timepoints
prop.test(Tol_contrasts.agree.Chl.ColTP$Agree, Tol_contrasts.agree.Chl.ColTP$Total)

Warning: Chi-squared approximation may be incorrect

    3-sample test for equality of proportions without continuity correction

data:  Tol_contrasts.agree.Chl.ColTP$Agree out of Tol_contrasts.agree.Chl.ColTP$Total
X-squared = 9, df = 2, p-value = 0.01111
alternative hypothesis: two.sided
sample estimates:
   prop 1    prop 2    prop 3 
0.8888889 1.0000000 0.4444444 

Sym vs Color Full

##Subset by Comparison
Tol_contrasts.agree.Sym.ColF<-subset(Tol_contrasts.agree.sum.tp, Compare=="Sym.ColF")

##Compare Proportions between Timepoints
prop.test(Tol_contrasts.agree.Sym.ColF$Agree, Tol_contrasts.agree.Sym.ColF$Total)

Warning: Chi-squared approximation may be incorrect

    3-sample test for equality of proportions without continuity correction

data:  Tol_contrasts.agree.Sym.ColF$Agree out of Tol_contrasts.agree.Sym.ColF$Total
X-squared = 3, df = 2, p-value = 0.2231
alternative hypothesis: two.sided
sample estimates:
   prop 1    prop 2    prop 3 
0.5555556 0.8888889 0.5555556 

Sym vs Color Timepoint

##Subset by Comparison
Tol_contrasts.agree.Sym.ColTP<-subset(Tol_contrasts.agree.sum.tp, Compare=="Sym.ColTP")

##Compare Proportions between Timepoints
prop.test(Tol_contrasts.agree.Sym.ColTP$Agree, Tol_contrasts.agree.Sym.ColTP$Total)

Warning: Chi-squared approximation may be incorrect

    3-sample test for equality of proportions without continuity correction

data:  Tol_contrasts.agree.Sym.ColTP$Agree out of Tol_contrasts.agree.Sym.ColTP$Total
X-squared = 5.0143, df = 2, p-value = 0.0815
alternative hypothesis: two.sided
sample estimates:
   prop 1    prop 2    prop 3 
0.6666667 1.0000000 0.5555556 

Color Full vs Color Timepoint

##Subset by Comparison
Tol_contrasts.agree.ColF.ColTP<-subset(Tol_contrasts.agree.sum.tp, Compare=="ColF.ColTP")

##Compare Proportions between Timepoints
prop.test(Tol_contrasts.agree.ColF.ColTP$Agree, Tol_contrasts.agree.ColF.ColTP$Total)

Warning: Chi-squared approximation may be incorrect

    3-sample test for equality of proportions without continuity correction

data:  Tol_contrasts.agree.ColF.ColTP$Agree out of Tol_contrasts.agree.ColF.ColTP$Total
X-squared = 1.08, df = 2, p-value = 0.5827
alternative hypothesis: two.sided
sample estimates:
   prop 1    prop 2    prop 3 
0.8888889 0.8888889 1.0000000 

Figures

Figure S2 Tolerance Across Metrics

Update for Panel

Tol_Chl_W2_SG.plot.fig<-ggplot(Tol_Chl_W2_SG, aes(x=Site.Geno, y=Chl.prop, colour=Genotype)) + 
  geom_errorbar(aes(ymin=Chl.prop-se, ymax=Chl.prop+se), width=cap.sz, linewidth=bar.sz)+
  geom_point(size=point.sz)+ 
  ggtitle("Summer 1")+
  theme_classic()+
  theme( axis.title.x = element_text(size = axis.title.sz), 
         axis.title.y = element_text(size = axis.title.sz), 
         axis.text.x=element_text(size=axis.txt.sz, colour="black"),
         axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
         legend.text=element_text(size=leg.txt.sz),
         legend.title=element_text(size=leg.title.sz), 
         legend.box.background = element_rect(color = "black"),
         legend.position="top", 
         legend.direction="horizontal",
         plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
  labs(x="", y="Proportion Chlorophyll Retained")+
  ylim(0, 1)+ 
  scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
  scale_color_manual(values = Geno.colors.o)+ 
  stat_pvalue_manual(data=Tol_Chl_W2_lm.p,  y.position=0.4, step.increase=0.25, label="Sig", hide.ns=TRUE)


Tol_Chl_M1_SG.plot.fig<-ggplot(Tol_Chl_M1_SG, aes(x=Site.Geno, y=Chl.prop, colour=Genotype)) + 
  geom_errorbar(aes(ymin=Chl.prop-se, ymax=Chl.prop+se), width=cap.sz, linewidth=bar.sz)+
  geom_point(size=point.sz)+ 
  ggtitle("Summer 2")+
  theme_classic()+
  theme( axis.title.x = element_text(size = axis.title.sz), 
         axis.title.y = element_text(size = axis.title.sz), 
         axis.text.x=element_text(size=axis.txt.sz, colour="black"),
         axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
         legend.text=element_text(size=leg.txt.sz),
         legend.title=element_text(size=leg.title.sz), 
         legend.box.background = element_rect(color = "black"),
         legend.position="top", 
         legend.direction="horizontal",
         plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
  labs(x="", y="")+
  ylim(0, 1)+ 
  scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
  scale_color_manual(values = Geno.colors.o)+ 
  stat_pvalue_manual(data=Tol_Chl_M1_lm.p,  y.position=0.4, step.increase=0.25, label="Sig", hide.ns=TRUE)


Tol_Chl_M4_SG.plot.fig<-ggplot(Tol_Chl_M4_SG, aes(x=Site.Geno, y=Chl.prop, colour=Genotype)) + 
  geom_errorbar(aes(ymin=Chl.prop-se, ymax=Chl.prop+se), width=cap.sz, linewidth=bar.sz)+
  geom_point(size=point.sz)+ 
  ggtitle("Winter")+
  theme_classic()+
  theme( axis.title.x = element_text(size = axis.title.sz), 
         axis.title.y = element_text(size = axis.title.sz), 
         axis.text.x=element_text(size=axis.txt.sz, colour="black"),
         axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
         legend.text=element_text(size=leg.txt.sz),
         legend.title=element_text(size=leg.title.sz), 
         legend.box.background = element_rect(color = "black"),
         legend.position="top", 
         legend.direction="horizontal",
         plot.title = element_text(size = plot.title.sz, colour="black", hjust = 0.5))+
  labs(x="", y="")+
  ylim(0, 1)+ 
  scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
  scale_color_manual(values = Geno.colors.o)+ 
  stat_pvalue_manual(data=Tol_Chl_M4_lm.p,  y.position=0.55, step.increase=0.20, label="Sig", hide.ns=TRUE)


Tol_Sym_W2_SG.plot.fig<-ggplot(Tol_Sym_W2_SG, aes(x=Site.Geno, y=Sym.prop, colour=Genotype)) + 
  geom_errorbar(aes(ymin=Sym.prop-se, ymax=Sym.prop+se), width=cap.sz, linewidth=bar.sz)+
  geom_point(size=point.sz)+ 
  theme_classic()+
  theme( axis.title.x = element_text(size = axis.title.sz), 
         axis.title.y = element_text(size = axis.title.sz), 
         axis.text.x=element_text(size=axis.txt.sz, colour="black"),
         axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
         legend.position="none")+
  labs(x="Site and Genotype", y="Proportion Symbionts Retained")+
  ylim(0, 1.5)+ 
  scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
  scale_color_manual(values = Geno.colors.o)+ 
  stat_pvalue_manual(data=Tol_Sym_W2_lm.p,  y.position=0.95, step.increase=0.15, label="Sig", hide.ns=TRUE)


Tol_Sym_M1_SG.plot.fig<-ggplot(Tol_Sym_M1_SG, aes(x=Site.Geno, y=Sym.prop, colour=Genotype)) + 
  geom_errorbar(aes(ymin=Sym.prop-se, ymax=Sym.prop+se), width=cap.sz, linewidth=bar.sz)+
  geom_point(size=point.sz)+ 
  theme_classic()+
  theme( axis.title.x = element_text(size = axis.title.sz), 
         axis.title.y = element_text(size = axis.title.sz), 
         axis.text.x=element_text(size=axis.txt.sz, colour="black"),
         axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
         legend.position="none")+
  labs(x="Site and Genotype", y="")+
  ylim(0, 1.5)+ 
  scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
  scale_color_manual(values = Geno.colors.o)+ 
  stat_pvalue_manual(data=Tol_Sym_M1_lm.p,  y.position=0.95, step.increase=0.15, label="Sig", hide.ns=TRUE)


Tol_Sym_M4_SG.plot.fig<-ggplot(Tol_Sym_M4_SG, aes(x=Site.Geno, y=Sym.prop, colour=Genotype)) + 
  geom_errorbar(aes(ymin=Sym.prop-se, ymax=Sym.prop+se), width=cap.sz, linewidth=bar.sz)+
  geom_point(size=point.sz)+ 
  theme_classic()+
  theme( axis.title.x = element_text(size = axis.title.sz), 
         axis.title.y = element_text(size = axis.title.sz), 
         axis.text.x=element_text(size=axis.txt.sz, colour="black"),
         axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
         legend.position="none")+
  labs(x="Site and Genotype", y="")+
  ylim(0, 1.5)+ 
  scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
  scale_color_manual(values = Geno.colors.o)+ 
  stat_pvalue_manual(data=Tol_Sym_M4_lm.p,  y.position=1.1, step.increase=0.15, label="Sig", hide.ns=TRUE)


Tol_ScoreF_W2_SG.plot.fig<-ggplot(Tol_ScoreF_W2_SG, aes(x=Site.Geno, y=Score_Full.prop, colour=Genotype)) + 
  geom_errorbar(aes(ymin=Score_Full.prop-se, ymax=Score_Full.prop+se), width=cap.sz, linewidth=bar.sz)+
  geom_point(size=point.sz)+ 
  theme_classic()+
  theme( axis.title.x = element_text(size = axis.title.sz), 
         axis.title.y = element_text(size = axis.title.sz), 
         axis.text.x=element_text(size=axis.txt.sz, colour="black"),
         axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
         legend.position="none")+
  labs(x="", y="Proportion Color Retained (Full)")+
  ylim(0.25, 1.5)+ 
  scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
  scale_color_manual(values = Geno.colors.o)+ 
  stat_pvalue_manual(data=Tol_ScoreF_W2_lm.p,  y.position=0.9, step.increase=0.35, label="Sig", hide.ns=TRUE)


Tol_ScoreF_M1_SG.plot.fig<-ggplot(Tol_ScoreF_M1_SG, aes(x=Site.Geno, y=Score_Full.prop, colour=Genotype)) + 
  geom_errorbar(aes(ymin=Score_Full.prop-se, ymax=Score_Full.prop+se), width=cap.sz, linewidth=bar.sz)+
  geom_point(size=point.sz)+ 
  theme_classic()+
  theme( axis.title.x = element_text(size = axis.title.sz), 
         axis.title.y = element_text(size = axis.title.sz), 
         axis.text.x=element_text(size=axis.txt.sz, colour="black"),
         axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
         legend.position="none")+
  labs(x="", y="")+
  ylim(0.25, 1.5)+ 
  scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
  scale_color_manual(values = Geno.colors.o)+ 
  stat_pvalue_manual(data=Tol_ScoreF_M1_lm.p,  y.position=1.1, step.increase=0.2, label="Sig", hide.ns=TRUE)


Tol_ScoreF_M4_SG.plot.fig<-ggplot(Tol_ScoreF_M4_SG, aes(x=Site.Geno, y=Score_Full.prop, colour=Genotype)) + 
  geom_errorbar(aes(ymin=Score_Full.prop-se, ymax=Score_Full.prop+se), width=cap.sz, linewidth=bar.sz)+
  geom_point(size=point.sz)+ 
  theme_classic()+
  theme( axis.title.x = element_text(size = axis.title.sz), 
         axis.title.y = element_text(size = axis.title.sz), 
         axis.text.x=element_text(size=axis.txt.sz, colour="black"),
         axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
         legend.position="none")+
  labs(x="", y="")+
  ylim(0.25, 1.5)+ 
  scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
  scale_color_manual(values = Geno.colors.o)+ 
  stat_pvalue_manual(data=Tol_ScoreF_M4_lm.p,  y.position=1.1, step.increase=0.22, label="Sig", hide.ns=TRUE)


Tol_ScoreTP_W2_SG.plot.fig<-ggplot(Tol_ScoreTP_W2_SG, aes(x=Site.Geno, y=Score_TP.prop, colour=Genotype)) + 
  geom_errorbar(aes(ymin=Score_TP.prop-se, ymax=Score_TP.prop+se), width=cap.sz, linewidth=bar.sz)+
  geom_point(size=point.sz)+ 
  theme_classic()+
  theme( axis.title.x = element_text(size = axis.title.sz), 
         axis.title.y = element_text(size = axis.title.sz), 
         axis.text.x=element_text(size=axis.txt.sz, colour="black"),
         axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
         legend.position="none")+
  labs(x="", y="Proportion Color Retained (Split)")+
  ylim(0.25, 1.5)+ 
  scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
  scale_color_manual(values = Geno.colors.o)+ 
  stat_pvalue_manual(data=Tol_ScoreTP_W2_lm.p,  y.position=1, step.increase=0.45, label="Sig", hide.ns=TRUE)


Tol_ScoreTP_M1_SG.plot.fig<-ggplot(Tol_ScoreTP_M1_SG, aes(x=Site.Geno, y=Score_TP.prop, colour=Genotype)) + 
  geom_errorbar(aes(ymin=Score_TP.prop-se, ymax=Score_TP.prop+se), width=cap.sz, linewidth=bar.sz)+
  geom_point(size=point.sz)+ 
  theme_classic()+
  theme( axis.title.x = element_text(size = axis.title.sz), 
         axis.title.y = element_text(size = axis.title.sz), 
         axis.text.x=element_text(size=axis.txt.sz, colour="black"),
         axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
         legend.position="none")+
  labs(x="", y="")+
  ylim(0.25, 1.5)+ 
  scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
  scale_color_manual(values = Geno.colors.o)+ 
  stat_pvalue_manual(data=Tol_ScoreTP_M1_lm.p,  y.position=1.05, step.increase=0.3, label="Sig", hide.ns=TRUE)


Tol_ScoreTP_M4_SG.plot.fig<-ggplot(Tol_ScoreTP_M4_SG, aes(x=Site.Geno, y=Score_TP.prop, colour=Genotype)) + 
  geom_errorbar(aes(ymin=Score_TP.prop-se, ymax=Score_TP.prop+se), width=cap.sz, linewidth=bar.sz)+
  geom_point(size=point.sz)+ 
  theme_classic()+
  theme( axis.title.x = element_text(size = axis.title.sz), 
         axis.title.y = element_text(size = axis.title.sz), 
         axis.text.x=element_text(size=axis.txt.sz, colour="black"),
         axis.text.y=element_text(size=axis.txt.sz, colour="black"), 
         legend.position="none")+
  labs(x="", y="")+
  ylim(0.25, 1.5)+ 
  scale_x_discrete(labels=c("","Klein","","","Something Special",""))+
  scale_color_manual(values = Geno.colors.o)+ 
  stat_pvalue_manual(data=Tol_ScoreTP_M4_lm.p,  y.position=1.1, step.increase=0.4, label="Sig", hide.ns=TRUE)

Figure S2

##Create Panel
Tolerance_fig<-plot_grid(Tol_Chl_W2_SG.plot.fig, Tol_Chl_M1_SG.plot.fig, Tol_Chl_M4_SG.plot.fig,
                  Tol_ScoreF_W2_SG.plot.fig, Tol_ScoreF_M1_SG.plot.fig, Tol_ScoreF_M4_SG.plot.fig,
                  Tol_ScoreTP_W2_SG.plot.fig, Tol_ScoreTP_M1_SG.plot.fig, Tol_ScoreTP_M4_SG.plot.fig,
                  Tol_Sym_W2_SG.plot.fig, Tol_Sym_M1_SG.plot.fig, Tol_Sym_M4_SG.plot.fig,
                        nrow=4, ncol=3, byrow=TRUE,
                        rel_widths=1,
                        rel_heights=c(1, 0.85, 0.85, 0.85), 
                        labels=c("A", "B", "C", 
                                 "D", "E", "F", 
                                 "G", "H", "I", 
                                 "J", "K", "L"), 
                        label_size=panel.lab.sz, 
                        label_fontface = "bold")

##Save Figure
ggsave(filename="Figures/FigureS2_Tolerance_Across_Metrics.png", plot=Tolerance_fig, dpi=300, width=14, height=20, units="in")

Figure 4 Contrasts Heatmap

##Save Figure
ggsave(filename="Figures/Figure4_Tolerance_Heatmap.png", plot=Tol_Pairs.plot, dpi=300, width=8, height=12, units="in")

Tables

Table S2 Linear Model Results

##Combine Results Tables
TableS2_Tol.lm.res<-rbind(Tol_Chl_W2_lm.res, Tol_Chl_M1_lm.res, Tol_Chl_M4_lm.res,
                        Tol_Sym_W2_lm.res, Tol_Sym_M1_lm.res, Tol_Sym_M4_lm.res,
                        Tol_ScoreF_W2_lm.res, Tol_ScoreF_M1_lm.res, Tol_ScoreF_M4_lm.res,
                    Tol_ScoreTP_W2_lm.res, Tol_ScoreTP_M1_lm.res, Tol_ScoreTP_M4_lm.res)

##Add a Season Variable
TableS2_Tol.lm.res$Season<-ifelse(TableS2_Tol.lm.res$TimeP=="W2", "Summer1", ifelse(TableS2_Tol.lm.res$TimeP=="M1", "Summer2", ifelse(TableS2_Tol.lm.res$TimeP=="M4", "Winter" , NA)))

##Organize
names(TableS2_Tol.lm.res)
TableS2_Tol.lm.res<-TableS2_Tol.lm.res[,c("TimeP", "Season", "Response", "Predictor", "DF", "Sum.Sq", "Mean.Sq", "F.value", "p.value", "EtaSq")]

#Round to 4 digits
TableS2_Tol.lm.res$Sum.Sq<-round(TableS2_Tol.lm.res$Sum.Sq, 4)
TableS2_Tol.lm.res$Mean.Sq<-round(TableS2_Tol.lm.res$Mean.Sq, 4)
TableS2_Tol.lm.res$F.value<-round(TableS2_Tol.lm.res$F.value, 4)
TableS2_Tol.lm.res$p.value<-round(TableS2_Tol.lm.res$p.value, 4)
TableS2_Tol.lm.res$EtaSq<-round(TableS2_Tol.lm.res$EtaSq, 4)

##Write Out Table
write.csv(TableS2_Tol.lm.res, "Tables/TableS2_Tolerance_LM_Results.csv", row.names=FALSE)

Table S3 Pairwise Comparison Results

##Pairwise Results Table
TableS3_Tol.pairwise<-Tol_contrasts

##Add a Season Variable
TableS3_Tol.pairwise$Season<-ifelse(TableS3_Tol.pairwise$TimeP=="W2", "Summer1", ifelse(TableS3_Tol.pairwise$TimeP=="M1", "Summer2", ifelse(TableS3_Tol.pairwise$TimeP=="M4", "Winter" , NA)))

##Re-order by Seasonal Timepoint
TableS3_Tol.pairwise<- TableS3_Tol.pairwise %>% arrange(as.numeric(TimeP))

##Organize
names(TableS3_Tol.pairwise)
TableS3_Tol.pairwise<-TableS3_Tol.pairwise[,c("TimeP", "Season", "Response", "Contrast", "Estimate", "SE", "df", "t.ratio", "p")]
TableS3_Tol.pairwise<-TableS3_Tol.pairwise %>% dplyr::rename( DF = df) %>% dplyr::rename( p.value = p)

#Round to 4 digits
TableS3_Tol.pairwise$Estimate<-round(TableS3_Tol.pairwise$Estimate, 4)
TableS3_Tol.pairwise$SE<-round(TableS3_Tol.pairwise$SE, 4)
TableS3_Tol.pairwise$t.ratio<-round(TableS3_Tol.pairwise$t.ratio, 4)
TableS3_Tol.pairwise$p.value<-round(TableS3_Tol.pairwise$p.value, 4)

##Write Out Table
write.csv(TableS3_Tol.pairwise, "Tables/TableS3_Tolerance_Pairwise_Results.csv", row.names=FALSE)